Code Packages:
# Load packages
library(gridExtra) # For using grid.table to save tables as images
## Warning: package 'gridExtra' was built under R version 4.4.3
library(flextable) # Alternative method to save df as image
## Warning: package 'flextable' was built under R version 4.4.3
library(tidyverse)
## Warning: package 'ggplot2' was built under R version 4.4.3
## Warning: package 'lubridate' was built under R version 4.4.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::combine() masks gridExtra::combine()
## ✖ purrr::compose() masks flextable::compose()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(zoo)
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(lme4)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
library(lmerTest)
##
## Attaching package: 'lmerTest'
##
## The following object is masked from 'package:lme4':
##
## lmer
##
## The following object is masked from 'package:stats':
##
## step
library(rlang)
##
## Attaching package: 'rlang'
##
## The following objects are masked from 'package:purrr':
##
## %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
## flatten_raw, invoke, splice
library(nortest) # For ad.test
library(car) # For homogeneity of variance with leveneTest
## Warning: package 'car' was built under R version 4.4.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.4.3
##
## Attaching package: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## some
library(rcompanion) # For scheirerRayHare test
## Warning: package 'rcompanion' was built under R version 4.4.3
library(pastecs) # For stat.desc
## Warning: package 'pastecs' was built under R version 4.4.3
##
## Attaching package: 'pastecs'
##
## The following objects are masked from 'package:dplyr':
##
## first, last
##
## The following object is masked from 'package:tidyr':
##
## extract
library(ggpubr) # For assumption plots
## Warning: package 'ggpubr' was built under R version 4.4.3
##
## Attaching package: 'ggpubr'
##
## The following objects are masked from 'package:flextable':
##
## border, font, rotate
library(conover.test) # For the conover.Iman() test function
library(DescTools) # Contain alternative ConoverTest function
## Warning: package 'DescTools' was built under R version 4.4.3
##
## Attaching package: 'DescTools'
##
## The following object is masked from 'package:car':
##
## Recode
library(rstatix) # tidyverse adapted stat tests and mshapiro_test() functions
## Warning: package 'rstatix' was built under R version 4.4.3
##
## Attaching package: 'rstatix'
##
## The following object is masked from 'package:stats':
##
## filter
library(DHARMa) # For multi-level lm models and shapiro test
## Warning: package 'DHARMa' was built under R version 4.4.3
## This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
library(patchwork) # For combining ggplots
## Warning: package 'patchwork' was built under R version 4.4.3
# Read and call data into df
# Primary samples df
ABL90 <- read.csv("Raw_Data/ABL90_Raw.csv")
# River's Parturition metadata
partuition_subcat <- read.csv("Raw_Data/River's_Rockfish_Metadata_Parturition_V7.csv")
# Data sifting: ABL90 dataset
# Step 1: Remove missing info
ABL_set1 <- ABL90 %>%
filter(Patient.ID_edited != "")
# Step 2: Separate Patient.ID into columns of sample types: blood plasma and instant freeze plasma
ABL_set2 <- separate(ABL_set1, 'Patient.ID_edited', into = c("Patient.ID_edited", "Sample_Type"), sep = ",")
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 4 rows [489, 619, 1034,
## 1177].
# Step 3:
# Convert Patient.ID_edited column data to character type
ABL_set2$Patient.ID_edited <- as.character(ABL_set2$Patient.ID_edited)
# Step 4 & 5: Sussing out specific sample errors
# Remove insufficient samples
ABL_set3 <- ABL_set2 %>%
filter(!is.na("Type")) %>%
filter(!str_detect(Errors.detected.during.measurement, "Insufficient sample"))
# Remove inhomogeneous samples
ABL_set3 <- ABL_set3 %>%
filter(!is.na("Type")) %>%
filter(!str_detect(Errors.detected.during.measurement, "Inhomogeneous sample"))
# Step 6: Filter for only blood samples
ABL_b_samp <- ABL_set3 %>%
filter(Sample_Type == "b")
Join ABL90 df with parturition_subcat df: Making ABL_merged
# Rename columns to match: Change 'ID' in parturition_subcat to 'Patient.ID_edited' so it is ready to join with ABL90 df
# Please note the 'Treatment' col in parturition_subcat does not match with the finalized 'Ambient.Or.OAH' col in metadata_atresia_guide or ABL90_merged df
# Rename 'ID' col to 'patient.ID_edited'
partuition_subcat <- partuition_subcat %>%
rename(Patient.ID_edited = ID)
# Connect main ABL90 dataset with my (River's) sub categorized parturition metadata
# ABL_merged <- partuition_subcat %>%
# left_join(ABL_b_samp, by = "Patient.ID_edited")
ABL_merged <- ABL_b_samp %>%
inner_join(partuition_subcat, by = "Patient.ID_edited")
Rename Columns Neatly:
# Rename Treatment & Parturition Outcome
ABL_merged <- ABL_merged %>%
rename(Ambient.Or.OAH = Consensus_General_Treatment,
Pregnant.Or.Atresia = Consensus_Atresia_Or_Pregnant)
Remove Replicate ID’s
# Cross Check all Columns: Cross validation to suss correct replicate row
check_params <- ABL_merged %>%
select(Patient.ID_edited, Time, Sample.., Measuring.Mode, Ambient.Or.OAH, Pregnant.Or.Atresia, pH, Glu..mmol.L., Hct...., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L., pCO2..mmHg.) %>%
filter(Patient.ID_edited == "9782D")
print(check_params)
## Patient.ID_edited Time Sample.. Measuring.Mode Ambient.Or.OAH
## 1 9782D 12/28/2023 10:04 1016 C 65uL Ambient
## 2 9782D 1/9/2024 17:21 863 S 65uL Ambient
## Pregnant.Or.Atresia pH Glu..mmol.L. Hct.... Na...mmol.L. Cl...mmol.L.
## 1 Post Parturition 7.348 0.5 9.8 184 149
## 2 Post Parturition 7.446 0.7 13.3 181 161
## K...mmol.L. Ca....mmol.L. pCO2..mmHg.
## 1 2.8 1.32 4.2
## 2 2.9 1.31 6.8
check_params <- ABL_merged %>%
select(Patient.ID_edited, Time, Sample.., Measuring.Mode, Ambient.Or.OAH, Pregnant.Or.Atresia, pH, Glu..mmol.L., Hct...., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L., pCO2..mmHg.) %>%
filter(Patient.ID_edited == "9783D")
print(check_params)
## Patient.ID_edited Time Sample.. Measuring.Mode Ambient.Or.OAH
## 1 9783D 1/19/2024 10:37 934 S 65uL Ambient
## 2 9783D 1/7/2024 12:32 860 S 65uL Ambient
## Pregnant.Or.Atresia pH Glu..mmol.L. Hct.... Na...mmol.L. Cl...mmol.L.
## 1 Post Parturition 7.484 1.7 0.2 175 193
## 2 Post Parturition 7.450 1.3 14.7 179 164
## K...mmol.L. Ca....mmol.L. pCO2..mmHg.
## 1 34.5 1.29 8.5
## 2 2.5 1.30 5.3
# Removed replicates: 9782D (2x) and 9783D (2x)
ABL_merged <- ABL_merged %>%
filter(Sample.. != "863",
Sample.. != "934")
Review Replicates:
| Row | ID | Time | Sample.. | Measuring.Mode | pH | Glu.mmol.L. | Hct…. | Na…mmol.L. | Cl…mmol.L. | K…mmol.L. | Ca.mmol.L. |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 32 | 9782D | 12/28/2023 10:04 | 1016 | C 65uL | 7.348 | 0.5 | 9.8 | 184 | 149 | 2.8 | 1.32 |
| 39 | 9782D | 1/9/2024 17:21 | 863 | S 65uL | 7.446 | 0.7 | 13.3 | 181 | 161 | 2.9 | 1.31 |
| 36 | 9783D | 1/19/2024 10:37 | 934 | S 65uL | 7.484 | 1.7 | 0.2 | 175 | 193 | 34.5 | 1.29 |
| 41 | 9783D | 1/7/2024 12:32 | 860 | S 65uL | 7.450 | 1.3 | 14.7 | 179 | 164 | 2.5 | 1.30 |
Methods for replicate sample removal: Check processing date of IDs, and remove samples that do not match that date.
Case ID_9782D:
Origianl ABL90 set contains:
9782D,b (Sample.. 1016, Time = 12/28/2023 10:04:00 AM)
9782D,p (Sample.. 1020, Time = 12/28/2023 10:25:00 AM)
9782D,p (Sample.. 865, Time = 1/9/2024 5:29:00 PM)
9782D,b (Sample.. 863, Time = 1/9/2024 5:21:00 PM)
2024_Fish_Metadata.csv says ID_9782D was processed 3/29/2024
Note Processing date in metadata will not match time date stamp in ABL90 data due to the machine being reset, and innacurrate dates being automatically recorded.
Deciding info: - For ID-9782D, refer to metadata sheet and look for other samples with that same processing date, then view thoes samples ABL90 ‘Time’. Try to find a matching ‘Time’ with other IDs and 9782D that were processed together during the same time recorded by the machine.
ID_97838,b (single obs, Sample.. 1024) Time = 12/28/2023 11:49:00 AM
Since 9782D,b Sample.. 1016 and 97838,b Sample.. 1024 both have a date timestamp occurrence on 12/28/2023, I have decided to keep 9782D,b Sample.. 1016 and drop 9783D,b Sample.. 863
Case ID_9783D
Original ABL90 dataset contains:
9783D,b (Sample.. 934, Time = 1/19/2024 10:37:00 AM)
9783D,c (Sample.. 861, Time = 1/7/2024 12:38:00 PM)
9783D,b (Sample.. 860, Time = 1/7/2024 12:32:00 PM)
2024_Fish_Metadata.csv says ID_9783D was processed 3/20/2024
Informed that ID_9783D: Sample.. 934 is a plasma sample (on account of its low hct value), remove this and keep Sample.. 860 (parameter values appear correct).
Remove Motalities and No info IDs and assign analysis ready data-frames:
# New df with Moralities removed: Note none of these samples made it into ABL90 df anyway, so looks like they are already filtered out.
Live_Samples <- ABL_merged %>%
filter(Patient.ID_edited != "9780C", # Mortality
Patient.ID_edited != "777AE", # Mortality
Patient.ID_edited != "777CA") # Mortality after parturition
# New df with mortality and 'No info' Id's removed
Primary_Samples <- Live_Samples %>%
filter(Patient.ID_edited != "777A0", # No info
Patient.ID_edited != "9782F", # No info
Patient.ID_edited != "777B3", # No info
Patient.ID_edited != "777AA", # No info
Patient.ID_edited != "777DE", # No info
Patient.ID_edited != "777CE", # No info
Patient.ID_edited != "777A6") # No info
# New df of Only Ambient Treatment: For testing parturition success
Ambient_Samples <- Primary_Samples %>%
filter(Ambient.Or.OAH == "Ambient")
# New df of Fecundity Samples:
Fecundity_Samples <- Primary_Samples %>%
filter(Fecundity_Or_Brood_Count != "NA",
Fecundity_Class != "NA") # removes 97706
# Fecundity samples without atresia Ids
Fecundity_No_Atresia_Samples <- Primary_Samples %>%
filter(All_Fecundity != "Na",
All_Fecundity != 0) # removes atresia IDs
# Fecundity Ambient Samples
Amb_Fec_Samples <- Primary_Samples %>%
filter(Ambient.Or.OAH == "Ambient",
All_Fecundity != "NA") # removes OAH treatment samples and ID 97706
# Fecundity Ambient Samples without atresia Ids
Amb_Fec_No_Atresia_Samples <- Primary_Samples %>%
filter(Ambient.Or.OAH == "Ambient", # removes OAH treatment
All_Fecundity != "NA", # removes all missing info IDs
All_Fecundity != 0) # removes atresia IDs
Change Data Types:
# Change Continuous Data to type to numeric, double, or factor
# Change data type of ions to factor or numeric
#glimpse(General_Samples)
#glimpse(Ambient_Samples)
Primary_Samples <- Primary_Samples %>%
mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))
Ambient_Samples <- Ambient_Samples %>%
mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))
Fecundity_Samples <- Fecundity_Samples %>%
mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))
Fecundity_No_Atresia_Samples <- Fecundity_Samples %>%
mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))
Amb_Fec_Samples <- Amb_Fec_Samples %>%
mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))
Amb_Fec_No_Atresia_Samples <- Amb_Fec_No_Atresia_Samples %>%
mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))
# glimpse(Primary_Samples)
# glimpse(Ambient_Samples)
# glimpse(Fecundity_Samples)
# glimpse(Amb_Fec_Samples)
# Fecundity data
#unique(Fecundity_Samples$Fecundity_Class)
# Arrange the order of parturition categories for visualizations & tests
# Primary Samples
# Arrange Treatment
Primary_Samples$Ambient.Or.OAH <- ordered(Primary_Samples$Ambient.Or.OAH, levels = c("Ambient", "OAH"))
# Arrange Parturition outcome
Primary_Samples$Pregnant.Or.Atresia <- ordered(Primary_Samples$Pregnant.Or.Atresia, levels = c("Post Parturition", "Atresia"))
# Arrange Brood Condition
Primary_Samples$Consensus_Brood_Condition <- ordered(Primary_Samples$Consensus_Brood_Condition, levels = c("Excellent", "Normal", "Low", "Very Low", "Atresia"))
# Arrange Fecundity Class
Primary_Samples$Fecundity_Class <- ordered(Primary_Samples$Fecundity_Class, levels = c("High (>50,000)", "Low (~1,000s)", "Atresia"))
# Remove NAs from ordered character factors
Primary_Samples <- Primary_Samples %>%
filter(Fecundity_Class != is.na(Fecundity_Class)) %>%
filter(Consensus_Brood_Condition != is.na(Consensus_Brood_Condition))
# Ambient data
Ambient_Samples$Ambient.Or.OAH <- ordered(Ambient_Samples$Ambient.Or.OAH, levels = c("Ambient", "OAH"))
Ambient_Samples$Pregnant.Or.Atresia <- ordered(Ambient_Samples$Pregnant.Or.Atresia, levels = c("Post Parturition", "Atresia"))
Ambient_Samples$Consensus_Brood_Condition <- ordered(Ambient_Samples$Consensus_Brood_Condition, levels = c("Excellent", "Normal", "Low", "Atresia"))
Ambient_Samples$Fecundity_Class <- ordered(Ambient_Samples$Fecundity_Class, levels = c("High (>50,000)", "Low (~1,000s)", "Atresia"))
# Remove NAs from ordered character factors
Ambient_Samples <- Ambient_Samples %>%
filter(Fecundity_Class != is.na(Fecundity_Class)) %>%
filter(Consensus_Brood_Condition != is.na(Consensus_Brood_Condition))
# Ambient Fecundity data
Amb_Fec_No_Atresia_Samples$Ambient.Or.OAH <- ordered(Amb_Fec_No_Atresia_Samples$Ambient.Or.OAH, levels = c("Ambient", "OAH"))
Amb_Fec_No_Atresia_Samples$Pregnant.Or.Atresia <- ordered(Amb_Fec_No_Atresia_Samples$Pregnant.Or.Atresia, levels = c("Post Parturition", "Atresia"))
Amb_Fec_No_Atresia_Samples$Consensus_Brood_Condition <- ordered(Amb_Fec_No_Atresia_Samples$Consensus_Brood_Condition, levels = c("Excellent", "Normal", "Low", "Atresia"))
Amb_Fec_No_Atresia_Samples$Fecundity_Class <- ordered(Amb_Fec_No_Atresia_Samples$Fecundity_Class, levels = c("High (>50,000)", "Low (~1,000s)", "Atresia"))
# Remove NAs from ordered character factors
Amb_Fec_No_Atresia_Samples <- Amb_Fec_No_Atresia_Samples %>%
filter(Fecundity_Class != is.na(Fecundity_Class)) %>%
filter(Consensus_Brood_Condition != is.na(Consensus_Brood_Condition))
# Save merged dataset in a worked folder
write.csv(x = Primary_Samples, file = "Worked-Data/Primary_Samples")
write.csv(x = Ambient_Samples, file = "Worked-Data/Ambient_Samples")
write.csv(x = Fecundity_Samples, file = "Worked-Data/Fecundity_Samples")
write.csv(x = Fecundity_No_Atresia_Samples, file = "Worked-Data/Fecundity_No_Atresia_Samples")
write.csv(x = Amb_Fec_Samples, file = "Worked-Data/Amb-Fec_Samples")
write.csv(x = Amb_Fec_No_Atresia_Samples, file = "Worked-Data/Amb_Fec_No_Atresia_Samples")
# Sample Size
# Ambient data:
ambient_sample_size_table <- Ambient_Samples %>%
group_by(Ambient.Or.OAH, Consensus_Brood_Condition) %>%
summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop")
print(ambient_sample_size_table)
## # A tibble: 4 × 3
## Ambient.Or.OAH Consensus_Brood_Condition n_size
## <ord> <ord> <int>
## 1 Ambient Excellent 1
## 2 Ambient Normal 4
## 3 Ambient Low 3
## 4 Ambient Atresia 13
pdf("data-images/ambient_sample_size_table.pdf")
grid.table(ambient_sample_size_table)
dev.off()
## png
## 2
png("data-images/ambient_sample_size_table.png")
grid.table(ambient_sample_size_table)
dev.off()
## png
## 2
ambient_fecundity_class_sample_size_table <- Ambient_Samples %>%
group_by(Ambient.Or.OAH, Fecundity_Class) %>%
summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop")
print(ambient_fecundity_class_sample_size_table)
## # A tibble: 3 × 3
## Ambient.Or.OAH Fecundity_Class n_size
## <ord> <ord> <int>
## 1 Ambient High (>50,000) 4
## 2 Ambient Low (~1,000s) 4
## 3 Ambient Atresia 13
pdf("data-images/ambient_fecundity_class_sample_size_table.pdf")
grid.table(ambient_fecundity_class_sample_size_table)
dev.off()
## png
## 2
png("data-images/ambient_fecundity_class_sample_size_table.png")
grid.table(ambient_fecundity_class_sample_size_table)
dev.off()
## png
## 2
Sample Size Barplot:
# View Sample Distributions
# Sample Size barplot
# Ambient Samples
# Consensus brood condition sample barplot
ambient_sample_size_barplot <- Ambient_Samples %>%
group_by(Consensus_Brood_Condition, Ambient.Or.OAH) %>%
summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop") %>%
ggplot(aes(x = Consensus_Brood_Condition, y = n_size)) +
geom_bar(stat = "identity") +
facet_grid(~ Ambient.Or.OAH) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 15.5), breaks = c(0, 3, 6, 9, 12, 15)) +
labs(title = "Ambient Sample Size",
x = "Parturition Outcome",
y = "Sample Size (n)") +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(ambient_sample_size_barplot)
ggsave(filename = "data-images/.png", plot = ambient_sample_size_barplot, device = "png")
## Saving 7 x 5 in image
ggsave(filename = "data-images/ambient_sample_size_barplot.pdf", plot = ambient_sample_size_barplot, device = "pdf")
## Saving 7 x 5 in image
# Fecundity class sample barplot
ambient_sample_size_barplot2 <- Ambient_Samples %>%
group_by(Fecundity_Class, Ambient.Or.OAH) %>%
summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop") %>%
ggplot(aes(x = Fecundity_Class, y = n_size)) +
geom_bar(stat = "identity") +
facet_grid(~ Ambient.Or.OAH) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 15.5), breaks = c(0, 3, 6, 9, 12, 15)) +
labs(title = "Ambient Sample Size",
x = "Parturition Outcome",
y = "Sample Size (n)") +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(ambient_sample_size_barplot2)
ggsave(filename = "data-images/.png", plot = ambient_sample_size_barplot2, device = "png")
## Saving 7 x 5 in image
ggsave(filename = "data-images/ambient_sample_size_barplot2.pdf", plot = ambient_sample_size_barplot2, device = "pdf")
## Saving 7 x 5 in image
# Amb Fec No atresia sample size
amb_fec_no_atresia_barplot <- Amb_Fec_No_Atresia_Samples %>%
group_by(Fecundity_Class, Ambient.Or.OAH) %>%
summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop") %>%
ggplot(aes(x = Fecundity_Class, y = n_size)) +
geom_bar(stat = "identity") +
facet_grid(~ Ambient.Or.OAH) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 15.5), breaks = c(4, 8, 12)) +
labs(title = "Ambient No Atresia Sample Size",
x = "Parturition Outcome",
y = "Sample Size (n)") +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(amb_fec_no_atresia_barplot)
ggsave(filename = "data-images/.png", plot = amb_fec_no_atresia_barplot, device = "png")
## Saving 7 x 5 in image
ggsave(filename = "data-images/amb_fec_no_atresia_barplot.pdf", plot = amb_fec_no_atresia_barplot, device = "pdf")
## Saving 7 x 5 in image
pH Summary Stats
# pH
# Summary stats
# Ambient data
Ambient_Samples %>%
group_by(Ambient.Or.OAH, Consensus_Brood_Condition) %>%
summarize(count = n(),
median = round(median(pH), 3),
mean = round(mean(pH), 3),
sd = round(sd(pH), 3),
cv = round(sd(pH)/mean(pH), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 4 × 7
## Ambient.Or.OAH Consensus_Brood_Condition count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Ambient Excellent 1 7.43 7.43 NA NA
## 2 Ambient Normal 4 7.48 7.45 0.087 0.012
## 3 Ambient Low 3 7.43 7.44 0.046 0.006
## 4 Ambient Atresia 13 7.37 7.38 0.1 0.014
Ambient_Samples %>%
group_by(Ambient.Or.OAH, Fecundity_Class) %>%
summarize(count = n(),
median = round(median(pH), 3),
mean = round(mean(pH), 3),
sd = round(sd(pH), 3),
cv = round(sd(pH)/mean(pH), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 3 × 7
## Ambient.Or.OAH Fecundity_Class count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Ambient High (>50,000) 4 7.44 7.43 0.076 0.01
## 2 Ambient Low (~1,000s) 4 7.46 7.45 0.054 0.007
## 3 Ambient Atresia 13 7.37 7.38 0.1 0.014
Amb_Fec_No_Atresia_Samples %>%
group_by(Ambient.Or.OAH, Fecundity_Class) %>%
summarize(count = n(),
median = round(median(pH), 3),
mean = round(mean(pH), 3),
sd = round(sd(pH), 3),
cv = round(sd(pH)/mean(pH), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 2 × 7
## Ambient.Or.OAH Fecundity_Class count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Ambient High (>50,000) 4 7.44 7.43 0.076 0.01
## 2 Ambient Low (~1,000s) 4 7.46 7.45 0.054 0.007
pH Boxplot
# pH boxplot
# Ambient Samples
pH_ambient_boxplot <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Consensus_Brood_Condition, y = pH)) +
geom_point(aes(x = Consensus_Brood_Condition, y = pH), alpha = 0.5, colour = "black") +
labs(title = "pH", x = "Parturition Success", y = "pH") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(pH_ambient_boxplot)
ggsave(filename = "data-images/pH_ambient_boxplot.pdf", plot = pH_ambient_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/pH_ambient_boxplot.png", plot = pH_ambient_boxplot, device = "png")
## Saving 7 x 5 in image
pH_ambient_fecundity_boxplot <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = pH)) +
geom_point(aes(x = Fecundity_Class, y = pH), alpha = 0.5, colour = "black") +
labs(title = "pH", x = "Fecundity Class", y = "pH") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(pH_ambient_fecundity_boxplot)
ggsave(filename = "data-images/pH_ambient_fecundity_boxplot.pdf", plot = pH_ambient_fecundity_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/pH_ambient_fecundity_boxplot.png", plot = pH_ambient_fecundity_boxplot, device = "png")
## Saving 7 x 5 in image
pH_amb_fec_no_atresia_boxplot <- ggplot(data = Amb_Fec_No_Atresia_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = pH)) +
geom_point(aes(x = Fecundity_Class, y = pH), alpha = 0.5, colour = "black") +
labs(title = "pH", x = "Fecundity Class", y = "pH") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(pH_amb_fec_no_atresia_boxplot)
ggsave(filename = "data-images/pH_amb_fec_no_atresia_boxplot.pdf", plot = pH_amb_fec_no_atresia_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/pH_amb_fec_no_atresia_boxplot.png", plot = pH_amb_fec_no_atresia_boxplot, device = "png")
## Saving 7 x 5 in image
pH_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
ggplot() +
geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = pH), colour = "black") +
labs(title = "pH", x = "Weight Adjusted Fecundity", y = "pH") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 20),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(pH_amb_fec_no_atresia_scatterplot)
ggsave(filename = "data-images/pH_amb_fec_no_atresia_scatterplot.pdf", plot = pH_amb_fec_no_atresia_scatterplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/pH_amb_fec_no_atresia_scatterplot.png", plot = pH_amb_fec_no_atresia_scatterplot, device = "png")
## Saving 7 x 5 in image
Differences: Between Brood Conditon
# pH
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
pH_aov_ambient_brood <- aov(pH ~ Consensus_Brood_Condition, data = Ambient_Samples)
summary(pH_aov_ambient_brood)
## Df Sum Sq Mean Sq F value Pr(>F)
## Consensus_Brood_Condition 3 0.02079 0.006930 0.795 0.513
## Residuals 17 0.14815 0.008715
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(pH_aov_ambient_brood)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = pH ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## $Consensus_Brood_Condition
## diff lwr upr p adj
## Normal-Excellent 0.01675000 -0.2799312 0.31343123 0.9984685
## Low-Excellent 0.00400000 -0.3024111 0.31041106 0.9999809
## Atresia-Excellent -0.05430769 -0.3296845 0.22106915 0.9423449
## Low-Normal -0.01275000 -0.2154219 0.18992187 0.9978869
## Atresia-Normal -0.07105769 -0.2227829 0.08066756 0.5566650
## Atresia-Low -0.05830769 -0.2282740 0.11165858 0.7651322
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$pH)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$pH
## W = 0.98534, p-value = 0.9808
shapiro.test(residuals(aov(pH ~ Consensus_Brood_Condition, data = Ambient_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov(pH ~ Consensus_Brood_Condition, data = Ambient_Samples))
## W = 0.95636, p-value = 0.446
# QQplot:
pH_ambient_res_qqplot <- ggqqplot(residuals(aov(pH ~ Consensus_Brood_Condition, data = Ambient_Samples))) +
labs(title = "Ambient pH Residual QQplot",
subtitle = "residuals(aov(pH ~ Consensus_Brood_Condition, data = Ambient_Samples))",
x = "pH Theoretical Quantiles (Predicted values)",
y = "pH Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(pH_ambient_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
#bartlett.test(pH ~ Consensus_Brood_Condition, data = Ambient_Samples)
# Nonparametric variance test (Levene's test):
leveneTest(pH ~ Consensus_Brood_Condition, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 0.671 0.5815
## 17
# Nonparametric stat test: Kruskall Wallis
kruskal.test(pH ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: pH by Consensus_Brood_Condition
## Kruskal-Wallis chi-squared = 2.8578, df = 3, p-value = 0.4141
# Nonparametric Post Hoc: Dunn's Test
DunnettTest(pH ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## Dunnett's test for comparing several treatments with a control :
## 95% family-wise confidence level
##
## $Excellent
## diff lwr.ci upr.ci pval
## Normal-Excellent 0.01675000 -0.2376275 0.2711275 0.9925
## Low-Excellent 0.00400000 -0.2587200 0.2667200 0.9999
## Atresia-Excellent -0.05430769 -0.2904186 0.1818032 0.8155
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pH Ambient Results
Which brood condition significantly different (if any)?
Parametric assumptions: - Normality (Shapiro-Wilk test): - Variance (Bartletts test): NA
Differences: Between Fecundity Class
# pH
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
pH_aov_ambient_fec <- aov(pH ~ Fecundity_Class, data = Ambient_Samples)
summary(pH_aov_ambient_fec)
## Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class 2 0.02183 0.010916 1.336 0.288
## Residuals 18 0.14711 0.008173
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(pH_aov_ambient_fec)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = pH ~ Fecundity_Class, data = Ambient_Samples)
##
## $Fecundity_Class
## diff lwr upr p adj
## Low (~1,000s)-High (>50,000) 0.02675000 -0.1363959 0.18989587 0.9084726
## Atresia-High (>50,000) -0.05080769 -0.1827287 0.08111329 0.5966404
## Atresia-Low (~1,000s) -0.07755769 -0.2094787 0.05436329 0.3141611
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$pH)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$pH
## W = 0.98534, p-value = 0.9808
shapiro.test(residuals(aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.96555, p-value = 0.6339
# QQplot:
pH_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(pH ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient pH Residual QQplot",
subtitle = "residuals(aov(pH ~ Fecundity_Class, data = Ambient_Samples))",
x = "pH Theoretical Quantiles (Predicted values)",
y = "pH Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(pH_ambient_fec_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(pH ~ Fecundity_Class, data = Ambient_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: pH by Fecundity_Class
## Bartlett's K-squared = 1.4521, df = 2, p-value = 0.4838
# Nonparametric variance test (Levene's test):
leveneTest(pH ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.44 0.6508
## 18
# Nonparametric stat test: Kruskall Wallis
kruskal.test(pH ~ Fecundity_Class, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: pH by Fecundity_Class
## Kruskal-Wallis chi-squared = 3.0375, df = 2, p-value = 0.219
# Nonparametric Post Hoc: Dunn's test
DunnettTest(pH ~ Fecundity_Class, data = Ambient_Samples)
##
## Dunnett's test for comparing several treatments with a control :
## 95% family-wise confidence level
##
## $`High (>50,000)`
## diff lwr.ci upr.ci pval
## Low (~1,000s)-High (>50,000) 0.02675000 -0.1252448 0.17874478 0.8717
## Atresia-High (>50,000) -0.05080769 -0.1737118 0.07209643 0.5071
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pH Ambient Results
Which fecundity class significantly different (if any)?
Parametric assumptions: - Normality (Shapiro-Wilk test): Met - Variance (Bartlett’s test): Met
Nonparametric variance test (Levene’s test):
Correlations: Parameter to Numerical Fecundity (without Atresia IDs)
# Ph
# Linear Model
pH_amb_fec_no_atresia_lm <- lm(pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples)
summary(pH_amb_fec_no_atresia_lm)
##
## Call:
## lm(formula = pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = Fecundity_No_Atresia_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.11940 -0.02163 0.01153 0.04429 0.05591
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 7.4752402 0.0369759 202.16
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight -0.0001977 0.0001939 -1.02
## Pr(>|t|)
## (Intercept) 4.01e-16 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.338
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.05852 on 8 degrees of freedom
## (18 observations deleted due to missingness)
## Multiple R-squared: 0.115, Adjusted R-squared: 0.004385
## F-statistic: 1.04 on 1 and 8 DF, p-value: 0.3377
pH_ambient_sim <- simulateResiduals(fittedModel = pH_amb_fec_no_atresia_lm)
plot(pH_ambient_sim)
## qu = 0.25, log(sigma) = -4.135132 : outer Newton did not converge fully.
plotResiduals(pH_amb_fec_no_atresia_lm)
## qu = 0.25, log(sigma) = -4.135132 : outer Newton did not converge fully.
testDispersion(pH_amb_fec_no_atresia_lm)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.90609, p-value = 0.952
## alternative hypothesis: two.sided
# Normality test
shapiro.test(residuals(pH_amb_fec_no_atresia_lm))
##
## Shapiro-Wilk normality test
##
## data: residuals(pH_amb_fec_no_atresia_lm)
## W = 0.89785, p-value = 0.2075
# Transform data: square root (undo log scale?)
pH_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
mutate(pH = sqrt(pH))
# Rerun lm model with square root transformed data
pH_amb_fec_no_atresia_lm_sqrt <- lm(pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = pH_ambient_sqrt)
summary(pH_amb_fec_no_atresia_lm_sqrt)
##
## Call:
## lm(formula = pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = pH_ambient_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.020889 -0.004920 0.002934 0.007526 0.010787
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 2.734e+00 7.412e-03 368.821
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight -4.078e-05 4.222e-05 -0.966
## Pr(>|t|)
## (Intercept) 2.68e-14 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.371
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01157 on 6 degrees of freedom
## Multiple R-squared: 0.1346, Adjusted R-squared: -0.009642
## F-statistic: 0.9331 on 1 and 6 DF, p-value: 0.3714
pH_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = pH_amb_fec_no_atresia_lm_sqrt)
plot(pH_amb_fec_sqrt_sim)
## qu = 0.25, log(sigma) = -2.642951 : outer Newton did not converge fully.
plotResiduals(pH_amb_fec_no_atresia_lm_sqrt)
## qu = 0.25, log(sigma) = -2.642951 : outer Newton did not converge fully.
testDispersion(pH_amb_fec_no_atresia_lm_sqrt)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
pH_amb_fec_no_atresia_lm_res <- residuals(pH_amb_fec_no_atresia_lm_sqrt)
# Normality test
shapiro.test(pH_amb_fec_no_atresia_lm_res) # normal
##
## Shapiro-Wilk normality test
##
## data: pH_amb_fec_no_atresia_lm_res
## W = 0.90414, p-value = 0.3147
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
#bartlett.test(pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples)
# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples)
Residuals Plot (if significant)
# pH
# Ambient samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
pH_ambient_lm_scatterplot <- ggplot(data = pH_ambient_sqrt, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = pH)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "pH", x = "Weight Adjusted Fecudnity", y = "Sqrt pH", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(pH_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/pH_ambient_lm_scatterplot.pdf", plot = pH_ambient_lm_scatterplot, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/pH_ambient_lm_scatterplot.png", plot = pH_ambient_lm_scatterplot, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
# sodium
# lm Boxplot
pH_ambient_lm_boxplot <- ggplot(data = pH_ambient_sqrt, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = pH)) +
geom_boxplot() +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "pH", x = "Weight Adjusted Fecundity", y = "Sqrt pH", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(pH_ambient_lm_boxplot)
## Warning: Continuous x aesthetic
## ℹ did you forget `aes(group = ...)`?
## `geom_smooth()` using formula = 'y ~ x'
hct Summary Stats
# hct
# Summary stats
# Ambient data
Ambient_Samples %>%
group_by(Ambient.Or.OAH, Consensus_Brood_Condition) %>%
summarize(count = n(),
median = round(median(Hct....), 3),
mean = round(mean(Hct....), 3),
sd = round(sd(Hct....), 3),
cv = round(sd(Hct....)/mean(Hct....), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 4 × 7
## Ambient.Or.OAH Consensus_Brood_Condition count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Ambient Excellent 1 18.2 18.2 NA NA
## 2 Ambient Normal 4 13.8 15.4 4.33 0.281
## 3 Ambient Low 3 14 13.4 3.29 0.245
## 4 Ambient Atresia 13 17.3 18.2 2.33 0.128
Ambient_Samples %>%
group_by(Ambient.Or.OAH, Fecundity_Class) %>%
summarize(count = n(),
median = round(median(Hct....), 3),
mean = round(mean(Hct....), 3),
sd = round(sd(Hct....), 3),
cv = round(sd(Hct....)/mean(Hct....), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 3 × 7
## Ambient.Or.OAH Fecundity_Class count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Ambient High (>50,000) 4 16.4 16.7 4.12 0.247
## 2 Ambient Low (~1,000s) 4 13.4 13.3 2.70 0.203
## 3 Ambient Atresia 13 17.3 18.2 2.33 0.128
Amb_Fec_No_Atresia_Samples %>%
group_by(Ambient.Or.OAH, Fecundity_Class) %>%
summarize(count = n(),
median = round(median(Hct....), 3),
mean = round(mean(Hct....), 3),
sd = round(sd(Hct....), 3),
cv = round(sd(Hct....)/mean(Hct....), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 2 × 7
## Ambient.Or.OAH Fecundity_Class count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Ambient High (>50,000) 4 16.4 16.7 4.12 0.247
## 2 Ambient Low (~1,000s) 4 13.4 13.3 2.70 0.203
hct Boxplot
# hct boxplot
# Ambient Samples
hct_ambient_boxplot <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Consensus_Brood_Condition, y = Hct....)) +
geom_point(aes(x = Consensus_Brood_Condition, y = Hct....), alpha = 0.5, colour = "black") +
labs(title = "Hematocrit", x = "Parturition Success", y = "Hematocrit (%)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(hct_ambient_boxplot)
ggsave(filename = "data-images/hct_ambient_boxplot.pdf", plot = hct_ambient_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/hct_ambient_boxplot.png", plot = hct_ambient_boxplot, device = "png")
## Saving 7 x 5 in image
hct_ambient_fecundity_boxplot <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Hct....)) +
geom_point(aes(x = Fecundity_Class, y = Hct....), alpha = 0.5, colour = "black") +
labs(title = "Hematocrit", x = "Fecundity Class", y = "Hematocrit (%)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(hct_ambient_fecundity_boxplot)
ggsave(filename = "data-images/hct_ambient_fecundity_boxplot.pdf", plot = hct_ambient_fecundity_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/hct_ambient_fecundity_boxplot.png", plot = hct_ambient_fecundity_boxplot, device = "png")
## Saving 7 x 5 in image
hct_amb_fec_no_atresia_boxplot <- ggplot(data = Amb_Fec_No_Atresia_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Hct....)) +
geom_point(aes(x = Fecundity_Class, y = Hct....), alpha = 0.5, colour = "black") +
labs(title = "Hematocrit", x = "Fecundity Class", y = "Hematocrit (%)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(hct_amb_fec_no_atresia_boxplot)
ggsave(filename = "data-images/hct_amb_fec_no_atresia_boxplot.pdf", plot = hct_amb_fec_no_atresia_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/hct_amb_fec_no_atresia_boxplot.png", plot = hct_amb_fec_no_atresia_boxplot, device = "png")
## Saving 7 x 5 in image
hct_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
ggplot() +
geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Hct....), colour = "black") +
labs(title = "Hematocrit", x = "Weight Adjusted Fecundity", y = "Hematocrit (%)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 20),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(hct_amb_fec_no_atresia_scatterplot)
ggsave(filename = "data-images/hct_amb_fec_no_atresia_scatterplot.pdf", plot = hct_amb_fec_no_atresia_scatterplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/hct_amb_fec_no_atresia_scatterplot.png", plot = hct_amb_fec_no_atresia_scatterplot, device = "png")
## Saving 7 x 5 in image
Differences: Between Brood Conditon
# hct
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
hct_aov_ambient_brood <- aov(Hct.... ~ Consensus_Brood_Condition, data = Ambient_Samples)
summary(hct_aov_ambient_brood) # Significant
## Df Sum Sq Mean Sq F value Pr(>F)
## Consensus_Brood_Condition 3 68.88 22.959 2.735 0.0758 .
## Residuals 17 142.69 8.394
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(hct_aov_ambient_brood) # Atresia-Low Significant
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Hct.... ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## $Consensus_Brood_Condition
## diff lwr upr p adj
## Normal-Excellent -2.825000e+00 -12.0325283 6.382528 0.8190431
## Low-Excellent -4.766667e+00 -14.2761610 4.742828 0.5017753
## Atresia-Excellent 3.552714e-15 -8.5463446 8.546345 1.0000000
## Low-Normal -1.941667e+00 -8.2316060 4.348273 0.8163755
## Atresia-Normal 2.825000e+00 -1.8838065 7.533807 0.3513279
## Atresia-Low 4.766667e+00 -0.5082517 10.041585 0.0843920
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$Hct....)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$Hct....
## W = 0.97482, p-value = 0.8355
shapiro.test(residuals(aov(Hct.... ~ Consensus_Brood_Condition, data = Ambient_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov(Hct.... ~ Consensus_Brood_Condition, data = Ambient_Samples))
## W = 0.9585, p-value = 0.4864
# QQplot:
hct_ambient_res_qqplot <- ggqqplot(residuals(aov(Hct.... ~ Consensus_Brood_Condition, data = Ambient_Samples))) +
labs(title = "Ambint Hematocrit Residual QQplot",
subtitle = "residuals(aov(Hct.... ~ Consensus_Brood_Condition, data = Ambient_Samples))",
x = "Hct Theoretical Quantiles (Predicted values)",
y = "Hct Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(hct_ambient_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
# bartlett.test(Hct.... ~ Consensus_Brood_Condition, data = Ambient_Samples)
# Nonparametric variance test (Levene's test):
leveneTest(Hct.... ~ Consensus_Brood_Condition, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 0.507 0.6826
## 17
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Hct.... ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Hct.... by Consensus_Brood_Condition
## Kruskal-Wallis chi-squared = 6.6608, df = 3, p-value = 0.08353
# Nonparametric Post Hoc: Dunn's test
DunnettTest(Hct.... ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## Dunnett's test for comparing several treatments with a control :
## 95% family-wise confidence level
##
## $Excellent
## diff lwr.ci upr.ci pval
## Normal-Excellent -2.825000 -10.719629 5.069629 0.5993
## Low-Excellent -4.766667 -12.920205 3.386871 0.2885
## Atresia-Excellent 0.000000 -7.327724 7.327724 1.0000
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Hematocrit Ambient Results
Which brood condition significantly different (if any)? - Atresia-Low (Tukey HSD: p adj = 0.0843920)
Parametric assumptions: - Normality (Shapiro-Wilk test): Met - Variance (Bartletts test): NA
Differences: Between Fecundity Class
# hct
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
hct_aov_ambient_fec <- aov(Hct.... ~ Fecundity_Class, data = Ambient_Samples)
summary(hct_aov_ambient_fec)
## Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class 2 73.83 36.92 4.824 0.021 *
## Residuals 18 137.74 7.65
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(hct_aov_ambient_fec)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Hct.... ~ Fecundity_Class, data = Ambient_Samples)
##
## $Fecundity_Class
## diff lwr upr p adj
## Low (~1,000s)-High (>50,000) -3.4 -8.3921454 1.592145 0.2186627
## Atresia-High (>50,000) 1.5 -2.5366864 5.536686 0.6176987
## Atresia-Low (~1,000s) 4.9 0.8633136 8.936686 0.0162654
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$Hct....)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$Hct....
## W = 0.97482, p-value = 0.8355
shapiro.test(residuals(aov(Hct.... ~ Fecundity_Class, data = Ambient_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov(Hct.... ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.97094, p-value = 0.7538
# QQplot:
hct_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(Hct.... ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Hematocrit Residual QQplot",
subtitle = "residuals(aov(Hct.... ~ Fecundity_Class, data = Ambient_Samples))",
x = "Hct Theoretical Quantiles (Predicted values)",
y = "Hct Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(hct_ambient_fec_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(Hct.... ~ Fecundity_Class, data = Ambient_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Hct.... by Fecundity_Class
## Bartlett's K-squared = 1.7219, df = 2, p-value = 0.4228
# Nonparametric variance test (Levene's test):
leveneTest(Hct.... ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 1.2003 0.3241
## 18
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Hct.... ~ Fecundity_Class, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Hct.... by Fecundity_Class
## Kruskal-Wallis chi-squared = 6.7821, df = 2, p-value = 0.03367
# Nonparametric Post Hoc: Dunn's Test
DunnettTest(Hct.... ~ Fecundity_Class, data = Ambient_Samples)
##
## Dunnett's test for comparing several treatments with a control :
## 95% family-wise confidence level
##
## $`High (>50,000)`
## diff lwr.ci upr.ci pval
## Low (~1,000s)-High (>50,000) -3.4 -8.050930 1.250930 0.1642
## Atresia-High (>50,000) 1.5 -2.260777 5.260777 0.5289
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Hct Ambient Results
Which fecundity class significantly different (if any)? - Atresia-Low (~1,000s): (Tukey HSD: p adj = 0.0162654)
Parametric assumptions: - Normality (Shapiro-Wilk test): Met - Variance (Bartlett’s test): Met
Nonparametric variance test (Levene’s test): Met
Correlations: Parameter to Numerical Fecundity (without Atresia IDs)
# hct
# Linear Model
hct_amb_fec_no_atresia_lm <- lm(Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples)
summary(hct_amb_fec_no_atresia_lm)
##
## Call:
## lm(formula = Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = Fecundity_No_Atresia_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.4692 -1.9674 -0.8183 1.4654 7.1131
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 13.09179 2.12136 6.171
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.01313 0.01112 1.181
## Pr(>|t|)
## (Intercept) 0.000268 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.271681
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.357 on 8 degrees of freedom
## (18 observations deleted due to missingness)
## Multiple R-squared: 0.1484, Adjusted R-squared: 0.04192
## F-statistic: 1.394 on 1 and 8 DF, p-value: 0.2717
hct_ambient_sim <- simulateResiduals(fittedModel = hct_amb_fec_no_atresia_lm)
plot(hct_ambient_sim)
## qu = 0.5, log(sigma) = -3.130892 : outer Newton did not converge fully.
## qu = 0.5, log(sigma) = -3.137393 : outer Newton did not converge fully.
## qu = 0.5, log(sigma) = -3.148809 : outer Newton did not converge fully.
## qu = 0.5, log(sigma) = -3.14559 : outer Newton did not converge fully.
plotResiduals(hct_amb_fec_no_atresia_lm)
## qu = 0.5, log(sigma) = -3.130892 : outer Newton did not converge fully.
## qu = 0.5, log(sigma) = -3.137393 : outer Newton did not converge fully.
## qu = 0.5, log(sigma) = -3.148809 : outer Newton did not converge fully.
## qu = 0.5, log(sigma) = -3.14559 : outer Newton did not converge fully.
testDispersion(hct_amb_fec_no_atresia_lm)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.90609, p-value = 0.952
## alternative hypothesis: two.sided
# Normality test
shapiro.test(residuals(hct_amb_fec_no_atresia_lm))
##
## Shapiro-Wilk normality test
##
## data: residuals(hct_amb_fec_no_atresia_lm)
## W = 0.89945, p-value = 0.2161
# Transform data: square root (undo log scale?)
hct_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
mutate(Hct.... = sqrt(Hct....))
# Rerun lm model with square root transformed data
hct_amb_fec_no_atresia_lm_sqrt <- lm(Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = hct_ambient_sqrt)
summary(hct_amb_fec_no_atresia_lm_sqrt)
##
## Call:
## lm(formula = Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = hct_ambient_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.45751 -0.18537 -0.07738 0.06130 0.87682
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 3.549914 0.292111 12.153
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.002034 0.001664 1.223
## Pr(>|t|)
## (Intercept) 1.89e-05 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.267
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.456 on 6 degrees of freedom
## Multiple R-squared: 0.1995, Adjusted R-squared: 0.06603
## F-statistic: 1.495 on 1 and 6 DF, p-value: 0.2673
hct_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = hct_amb_fec_no_atresia_lm_sqrt)
plot(hct_amb_fec_sqrt_sim)
plotResiduals(hct_amb_fec_no_atresia_lm_sqrt)
testDispersion(hct_amb_fec_no_atresia_lm_sqrt)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
hct_amb_fec_no_atresia_lm_res <- residuals(hct_amb_fec_no_atresia_lm_sqrt)
# Normality test
shapiro.test(hct_amb_fec_no_atresia_lm_res)
##
## Shapiro-Wilk normality test
##
## data: hct_amb_fec_no_atresia_lm_res
## W = 0.87306, p-value = 0.1615
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
#bartlett.test(Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = hct_ambient_sqrt)
# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = hct_ambient_sqrt)
Hct Ambient No Atresia Results
DHARMA Warning: Data dispersion error
Parametric assumptions: - Normality (Shapiro-Wilk test): Met - Variance (Bartlett’s test): NA
Nonparametric variance test (Levene’s test): NA
Residuals Plot (if significant)
# hct
# Ambient samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
hct_ambient_lm_scatterplot <- ggplot(data = hct_ambient_sqrt, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Hct....)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Hematocrit", x = "Weight Adjusted Fecudnity", y = "Sqrt Hematocrit (%)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(hct_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/hct_ambient_lm_scatterplot.pdf", plot = hct_ambient_lm_scatterplot, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
ggsave(filename = "data-images/hct_ambient_lm_scatterplot.png", plot = hct_ambient_lm_scatterplot, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
glucose Summary Stats
# glucose
# Summary stats
# Ambient data
Ambient_Samples %>%
group_by(Ambient.Or.OAH, Consensus_Brood_Condition) %>%
summarize(count = n(),
median = round(median(Glu..mmol.L.), 3),
mean = round(mean(Glu..mmol.L.), 3),
sd = round(sd(Glu..mmol.L.), 3),
cv = round(sd(Glu..mmol.L.)/mean(Glu..mmol.L.), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 4 × 7
## Ambient.Or.OAH Consensus_Brood_Condition count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Ambient Excellent 1 1.8 1.8 NA NA
## 2 Ambient Normal 4 1.55 1.58 0.275 0.175
## 3 Ambient Low 3 1.1 1.6 1.04 0.653
## 4 Ambient Atresia 13 2 2.1 0.698 0.332
Ambient_Samples %>%
group_by(Ambient.Or.OAH, Fecundity_Class) %>%
summarize(count = n(),
median = round(median(Glu..mmol.L.), 3),
mean = round(mean(Glu..mmol.L.), 3),
sd = round(sd(Glu..mmol.L.), 3),
cv = round(sd(Glu..mmol.L.)/mean(Glu..mmol.L.), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 3 × 7
## Ambient.Or.OAH Fecundity_Class count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Ambient High (>50,000) 4 1.55 1.55 0.238 0.154
## 2 Ambient Low (~1,000s) 4 1.5 1.68 0.866 0.517
## 3 Ambient Atresia 13 2 2.1 0.698 0.332
Amb_Fec_No_Atresia_Samples %>%
group_by(Ambient.Or.OAH, Fecundity_Class) %>%
summarize(count = n(),
median = round(median(Glu..mmol.L.), 3),
mean = round(mean(Glu..mmol.L.), 3),
sd = round(sd(Glu..mmol.L.), 3),
cv = round(sd(Glu..mmol.L.)/mean(Glu..mmol.L.), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 2 × 7
## Ambient.Or.OAH Fecundity_Class count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Ambient High (>50,000) 4 1.55 1.55 0.238 0.154
## 2 Ambient Low (~1,000s) 4 1.5 1.68 0.866 0.517
Glucose Boxplot
# glucose boxplot
# Ambient Samples
glucose_ambient_boxplot <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Consensus_Brood_Condition, y = Glu..mmol.L.)) +
geom_point(aes(x = Consensus_Brood_Condition, y = Glu..mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Glucose", x = "Parturition Success", y = "Glucose (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(glucose_ambient_boxplot)
ggsave(filename = "data-images/glucose_ambient_boxplot.pdf", plot = glucose_ambient_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/glucose_ambient_boxplot.png", plot = glucose_ambient_boxplot, device = "png")
## Saving 7 x 5 in image
glucose_ambient_fecundity_boxplot <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Glu..mmol.L.)) +
geom_point(aes(x = Fecundity_Class, y = Glu..mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Glucose", x = "Fecundity Class", y = "Glucose (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(glucose_ambient_fecundity_boxplot)
ggsave(filename = "data-images/glucose_ambient_fecundity_boxplot.pdf", plot = glucose_ambient_fecundity_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/glucose_ambient_fecundity_boxplot.png", plot = glucose_ambient_fecundity_boxplot, device = "png")
## Saving 7 x 5 in image
glucose_amb_fec_no_atresia_boxplot <- ggplot(data = Amb_Fec_No_Atresia_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Glu..mmol.L.)) +
geom_point(aes(x = Fecundity_Class, y = Glu..mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Glucose", x = "Fecundity Class", y = "Glucose (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(glucose_amb_fec_no_atresia_boxplot)
ggsave(filename = "data-images/glucose_amb_fec_no_atresia_boxplot.pdf", plot = glucose_amb_fec_no_atresia_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/glucose_amb_fec_no_atresia_boxplot.png", plot = glucose_amb_fec_no_atresia_boxplot, device = "png")
## Saving 7 x 5 in image
glucose_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
ggplot() +
geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Glu..mmol.L.), colour = "black") +
labs(title = "Glucose", x = "Weight Adjusted Fecundity", y = "Glucose (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 20),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(glucose_amb_fec_no_atresia_scatterplot)
ggsave(filename = "data-images/glucose_amb_fec_no_atresia_scatterplot.pdf", plot = glucose_amb_fec_no_atresia_scatterplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/glucose_amb_fec_no_atresia_scatterplot.pdf", plot = glucose_amb_fec_no_atresia_scatterplot, device = "png")
## Saving 7 x 5 in image
Differences: Between Brood Conditon
# glucose
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
glucose_aov_ambient_brood <- aov(Glu..mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
summary(glucose_aov_ambient_brood)
## Df Sum Sq Mean Sq F value Pr(>F)
## Consensus_Brood_Condition 3 1.218 0.4061 0.837 0.492
## Residuals 17 8.247 0.4851
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(glucose_aov_ambient_brood)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Glu..mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## $Consensus_Brood_Condition
## diff lwr upr p adj
## Normal-Excellent -0.225 -2.4386100 1.988610 0.9912859
## Low-Excellent -0.200 -2.4862066 2.086207 0.9943918
## Atresia-Excellent 0.300 -1.7546528 2.354653 0.9751438
## Low-Normal 0.025 -1.4871835 1.537184 0.9999611
## Atresia-Normal 0.525 -0.6070586 1.657059 0.5643929
## Atresia-Low 0.500 -0.7681593 1.768159 0.6821448
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$Glu..mmol.L.)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$Glu..mmol.L.
## W = 0.85634, p-value = 0.005473
shapiro.test(residuals(aov(Glu..mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov(Glu..mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))
## W = 0.79895, p-value = 0.0006323
# QQplot:
glucose_ambient_res_qqplot <- ggqqplot(residuals(aov(Glu..mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))) +
labs(title = "Glucose Interactive Residual QQplot",
subtitle = "residuals(aov(Glu..mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))",
x = "Glucose Theoretical Quantiles (Predicted values)",
y = "Glucose Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(glucose_ambient_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
# bartlett.test(Glu..mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
# Nonparametric variance test (Levene's test):
leveneTest(Glu..mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 0.4505 0.7202
## 17
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Glu..mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Glu..mmol.L. by Consensus_Brood_Condition
## Kruskal-Wallis chi-squared = 3.9713, df = 3, p-value = 0.2646
# Nonparametric Post Hoc: Dunn's test
DunnTest(Glu..mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## Normal-Excellent -3.5000000 1.0000
## Low-Excellent -2.8333333 1.0000
## Atresia-Excellent 2.5384615 1.0000
## Low-Normal 0.6666667 1.0000
## Atresia-Normal 6.0384615 0.5244
## Atresia-Low 5.3717949 0.8733
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Glucose Ambient Results
Parametric assumptions: - Normality (Shapiro-Wilk test): Not met - Variance (Bartletts test): NA
Differences: Between Fecundity Class
# glucose
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
glucose_aov_ambient_fec <- aov(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
summary(glucose_aov_ambient_fec)
## Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class 2 1.208 0.6041 1.317 0.293
## Residuals 18 8.257 0.4587
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(glucose_aov_ambient_fec)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## $Fecundity_Class
## diff lwr upr p adj
## Low (~1,000s)-High (>50,000) 0.125 -1.0973103 1.347310 0.9632210
## Atresia-High (>50,000) 0.550 -0.4383693 1.538369 0.3518867
## Atresia-Low (~1,000s) 0.425 -0.5633693 1.413369 0.5277964
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$Glu..mmol.L.)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$Glu..mmol.L.
## W = 0.85634, p-value = 0.005473
shapiro.test(residuals(aov(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.82056, p-value = 0.001376
# QQplot:
glucose_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Glucose Interactive Residual QQplot",
subtitle = "residuals(aov(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples))",
x = "Glucose Theoretical Quantiles (Predicted values)",
y = "Glucose Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(glucose_ambient_fec_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Glu..mmol.L. by Fecundity_Class
## Bartlett's K-squared = 3.669, df = 2, p-value = 0.1597
# Nonparametric variance test (Levene's test):
leveneTest(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.9391 0.4093
## 18
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Glu..mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 4.0755, df = 2, p-value = 0.1303
# Nonparametric Post Hoc: Dunn's test
DunnettTest(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Dunnett's test for comparing several treatments with a control :
## 95% family-wise confidence level
##
## $`High (>50,000)`
## diff lwr.ci upr.ci pval
## Low (~1,000s)-High (>50,000) 0.125 -1.0137649 1.263765 0.9470
## Atresia-High (>50,000) 0.550 -0.3708139 1.470814 0.2758
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Glucose Ambient Results
Parametric assumptions: - Normality (Shapiro-Wilk test): Not met - Variance (Bartlett’s test): Met
Nonparametric variance test (Levene’s test): Met
Correlations: Parameter to Numerical Fecundity (without Atresia IDs)
# glucose
# Linear Model
glucose_amb_fec_no_atresia_lm <- lm(Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples)
summary(glucose_amb_fec_no_atresia_lm)
##
## Call:
## lm(formula = Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = Fecundity_No_Atresia_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.8882 -0.5589 -0.3191 0.2610 1.7458
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.352797 0.563859 2.399
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.002951 0.002957 0.998
## Pr(>|t|)
## (Intercept) 0.0432 *
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.3475
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8924 on 8 degrees of freedom
## (18 observations deleted due to missingness)
## Multiple R-squared: 0.1107, Adjusted R-squared: -0.0004362
## F-statistic: 0.9961 on 1 and 8 DF, p-value: 0.3475
glucose_ambient_sim <- simulateResiduals(fittedModel = glucose_amb_fec_no_atresia_lm)
plot(glucose_ambient_sim)
plotResiduals(glucose_amb_fec_no_atresia_lm)
testDispersion(glucose_amb_fec_no_atresia_lm)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.90609, p-value = 0.952
## alternative hypothesis: two.sided
# Normality test:
shapiro.test(residuals(glucose_amb_fec_no_atresia_lm))
##
## Shapiro-Wilk normality test
##
## data: residuals(glucose_amb_fec_no_atresia_lm)
## W = 0.86864, p-value = 0.09638
# Transform data: square root
glucose_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
mutate(Glu..mmol.L. = sqrt(Glu..mmol.L.))
# Rerun lm model with square root transformed data
glucose_amb_fec_no_atresia_lm_sqrt <- lm(Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = glucose_ambient_sqrt)
summary(glucose_amb_fec_no_atresia_lm_sqrt)
##
## Call:
## lm(formula = Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = glucose_ambient_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.30380 -0.13647 -0.01130 0.07862 0.42723
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.2241905 0.1551624 7.890
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.0001917 0.0008838 0.217
## Pr(>|t|)
## (Intercept) 0.00022 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.83544
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2422 on 6 degrees of freedom
## Multiple R-squared: 0.007783, Adjusted R-squared: -0.1576
## F-statistic: 0.04706 on 1 and 6 DF, p-value: 0.8354
glucose_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = glucose_amb_fec_no_atresia_lm_sqrt)
plot(glucose_amb_fec_sqrt_sim)
plotResiduals(glucose_amb_fec_no_atresia_lm_sqrt)
testDispersion(glucose_amb_fec_no_atresia_lm_sqrt)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
glucose_amb_fec_no_atresia_lm_res <- residuals(glucose_amb_fec_no_atresia_lm_sqrt)
# Normality test
shapiro.test(glucose_amb_fec_no_atresia_lm_res)
##
## Shapiro-Wilk normality test
##
## data: glucose_amb_fec_no_atresia_lm_res
## W = 0.96004, p-value = 0.8104
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
#bartlett.test(Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = glucose_ambient_sqrt)
# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = glucose_ambient_sqrt)
Glucose Ambient No Atresia Results
DHARMA Warning: Data dispersion error
Parametric assumptions: - Normality (Shapiro-Wilk test): Met ish - Variance (Bartlett’s test): NA
Nonparametric variance test (Levene’s test): NA
Residuals Plot (if significant)
# glucose
# Ambient samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
glucose_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Glu..mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Glucose", x = "Weight Adjusted Fecudnity", y = "Glucose (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(glucose_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave(filename = "data-images/glucose_ambient_lm_scatterplot.pdf", plot = glucose_ambient_lm_scatterplot, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave(filename = "data-images/glucose_ambient_lm_scatterplot.png", plot = glucose_ambient_lm_scatterplot, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
Sodium Summary Stats
# sodium
# Summary stats
# Ambient data
Ambient_Samples %>%
group_by(Ambient.Or.OAH, Consensus_Brood_Condition) %>%
summarize(count = n(),
median = round(median(Na...mmol.L.), 3),
mean = round(mean(Na...mmol.L.), 3),
sd = round(sd(Na...mmol.L.), 3),
cv = round(sd(Na...mmol.L.)/mean(Na...mmol.L.), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 4 × 7
## Ambient.Or.OAH Consensus_Brood_Condition count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Ambient Excellent 1 177 177 NA NA
## 2 Ambient Normal 4 178. 177. 2.87 0.016
## 3 Ambient Low 3 175 173. 3.79 0.022
## 4 Ambient Atresia 13 178 181. 11.4 0.063
Ambient_Samples %>%
group_by(Ambient.Or.OAH, Fecundity_Class) %>%
summarize(count = n(),
median = round(median(Na...mmol.L.), 3),
mean = round(mean(Na...mmol.L.), 3),
sd = round(sd(Na...mmol.L.), 3),
cv = round(sd(Na...mmol.L.)/mean(Na...mmol.L.), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 3 × 7
## Ambient.Or.OAH Fecundity_Class count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Ambient High (>50,000) 4 178. 177. 2.63 0.015
## 2 Ambient Low (~1,000s) 4 176. 175. 4.19 0.024
## 3 Ambient Atresia 13 178 181. 11.4 0.063
Amb_Fec_No_Atresia_Samples %>%
group_by(Ambient.Or.OAH, Fecundity_Class) %>%
summarize(count = n(),
median = round(median(Na...mmol.L.), 3),
mean = round(mean(Na...mmol.L.), 3),
sd = round(sd(Na...mmol.L.), 3),
cv = round(sd(Na...mmol.L.)/mean(Na...mmol.L.), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 2 × 7
## Ambient.Or.OAH Fecundity_Class count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Ambient High (>50,000) 4 178. 177. 2.63 0.015
## 2 Ambient Low (~1,000s) 4 176. 175. 4.19 0.024
Sodium Boxplot
# Na+ boxplot
# Ambient Samples
sodium_ambient_boxplot <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Consensus_Brood_Condition, y = Na...mmol.L.)) +
geom_point(aes(x = Consensus_Brood_Condition, y = Na...mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Sodium", x = "Parturition Success", y = "Sodium (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(sodium_ambient_boxplot)
ggsave(filename = "data-images/sodium_ambient_boxplot.pdf", plot = sodium_ambient_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/sodium_ambient_boxplot.png", plot = sodium_ambient_boxplot, device = "png")
## Saving 7 x 5 in image
sodium_ambient_fecundity_boxplot <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Na...mmol.L.)) +
geom_point(aes(x = Fecundity_Class, y = Na...mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Sodium", x = "Fecundity Class", y = "Sodium (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(sodium_ambient_fecundity_boxplot)
ggsave(filename = "data-images/sodium_ambient_fecundity_boxplot.pdf", plot = sodium_ambient_fecundity_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/sodium_ambient_fecundity_boxplot.png", plot = sodium_ambient_fecundity_boxplot, device = "png")
## Saving 7 x 5 in image
sodium_amb_fec_no_atresia_boxplot <- ggplot(data = Amb_Fec_No_Atresia_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Na...mmol.L.)) +
geom_point(aes(x = Fecundity_Class, y = Na...mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Sodium", x = "Fecundity Class", y = "Sodium (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(sodium_amb_fec_no_atresia_boxplot)
ggsave(filename = "data-images/sodium_amb_fec_no_atresia_boxplot.pdf", plot = sodium_amb_fec_no_atresia_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/sodium_amb_fec_no_atresia_boxplot.png", plot = sodium_amb_fec_no_atresia_boxplot, device = "png")
## Saving 7 x 5 in image
sodium_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
ggplot() +
geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Na...mmol.L.), colour = "black") +
labs(title = "Sodium", x = "Weight Adjusted Fecundity", y = "Sodium (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 20),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(sodium_amb_fec_no_atresia_scatterplot)
ggsave(filename = "data-images/sodium_amb_fec_no_atresia_scatterplot.pdf", plot = sodium_amb_fec_no_atresia_scatterplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/sodium_amb_fec_no_atresia_scatterplot.pdf", plot = sodium_amb_fec_no_atresia_scatterplot, device = "png")
## Saving 7 x 5 in image
Differences: Between Brood Conditon
# Na+
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
sodium_aov_ambient_brood <- aov(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
summary(sodium_aov_ambient_brood)
## Df Sum Sq Mean Sq F value Pr(>F)
## Consensus_Brood_Condition 3 176.8 58.95 0.622 0.611
## Residuals 17 1611.7 94.81
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(sodium_aov_ambient_brood)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## $Consensus_Brood_Condition
## diff lwr upr p adj
## Normal-Excellent 0.250000 -30.694633 31.19463 0.9999955
## Low-Excellent -3.666667 -35.626146 28.29281 0.9875929
## Atresia-Excellent 4.230769 -24.491760 32.95330 0.9745138
## Low-Normal -3.916667 -25.055875 17.22254 0.9514394
## Atresia-Normal 3.980769 -11.844573 19.80611 0.8897586
## Atresia-Low 7.897436 -9.830494 25.62537 0.5954090
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$Na...mmol.L.)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$Na...mmol.L.
## W = 0.52676, p-value = 3.625e-07
shapiro.test(residuals(aov(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))
## W = 0.54475, p-value = 5.408e-07
# QQplot:
sodium_ambient_res_qqplot <- ggqqplot(residuals(aov(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))) +
labs(title = "Sodium Residual QQplot",
subtitle = "residuals(aov(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))",
x = "Sodium Theoretical Quantiles (Predicted values)",
y = "Sodium Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(sodium_ambient_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
# bartlett.test(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
# Nonparametric variance test (Levene's test):
leveneTest(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 0.194 0.899
## 17
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Na...mmol.L. by Consensus_Brood_Condition
## Kruskal-Wallis chi-squared = 5.5569, df = 3, p-value = 0.1353
# Nonparametric Post Hoc: Dunns
DunnTest(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## Normal-Excellent 3.0000000 1.0000
## Low-Excellent -5.5000000 1.0000
## Atresia-Excellent 3.5769231 1.0000
## Low-Normal -8.5000000 0.3481
## Atresia-Normal 0.5769231 1.0000
## Atresia-Low 9.0769231 0.1252
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sodium Ambient Results
Parametric assumptions: - Normality (Shapiro-Wilk test): Not met - Variance (Bartletts test): NA
Nonparametric variance (Levene’s test): Met
Nonparametric Post Hoc: Dunn’s Test
Atresia-Low: mean.rank.diff = 9.0769231, pval = 0.1252
Differences: Between Fecundity Class
# Na+
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
sodium_aov_ambient_fec <- aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
summary(sodium_aov_ambient_fec)
## Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class 2 156.8 78.38 0.865 0.438
## Residuals 18 1631.8 90.66
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(sodium_aov_ambient_fec)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## $Fecundity_Class
## diff lwr upr p adj
## Low (~1,000s)-High (>50,000) -2.000000 -19.18271 15.18271 0.9526465
## Atresia-High (>50,000) 4.480769 -9.41330 18.37484 0.6939179
## Atresia-Low (~1,000s) 6.480769 -7.41330 20.37484 0.4737484
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$Na...mmol.L.)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$Na...mmol.L.
## W = 0.52676, p-value = 3.625e-07
shapiro.test(residuals(aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.56321, p-value = 8.235e-07
# QQplot:
sodium_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Sodium Residual QQplot",
subtitle = "residuals(aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))",
x = "Sodium Theoretical Quantiles (Predicted values)",
y = "Sodium Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(sodium_ambient_fec_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Na...mmol.L. by Fecundity_Class
## Bartlett's K-squared = 7.4634, df = 2, p-value = 0.02395
# Nonparametric variance test (Levene's test):
leveneTest(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.2098 0.8127
## 18
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Na...mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 2.8608, df = 2, p-value = 0.2392
# Nonparametric Post Hoc: Dunn's test
DunnettTest(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Dunnett's test for comparing several treatments with a control :
## 95% family-wise confidence level
##
## $`High (>50,000)`
## diff lwr.ci upr.ci pval
## Low (~1,000s)-High (>50,000) -2.000000 -18.008265 14.00826 0.9321
## Atresia-High (>50,000) 4.480769 -8.463634 17.42517 0.6108
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sodium Ambient Results
Parametric assumptions: - Normality (Shapiro-Wilk test): Not met - Variance (Bartlett’s test): Not Met
Nonparametric variance test (Levene’s test): Met
Correlations: Parameter to Numerical Fecundity (without Atresia IDs)
# Na+
# Linear Model
sodium_amb_fec_no_atresia_lm <- lm(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples)
summary(sodium_amb_fec_no_atresia_lm)
##
## Call:
## lm(formula = Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = Fecundity_No_Atresia_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.1021 -0.3010 0.6449 1.9545 3.3374
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 175.27869 2.08086 84.234
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.00558 0.01091 0.511
## Pr(>|t|)
## (Intercept) 4.4e-13 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.623
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.293 on 8 degrees of freedom
## (18 observations deleted due to missingness)
## Multiple R-squared: 0.03166, Adjusted R-squared: -0.08939
## F-statistic: 0.2615 on 1 and 8 DF, p-value: 0.6229
sodium_ambient_sim <- simulateResiduals(fittedModel = sodium_amb_fec_no_atresia_lm)
plot(sodium_ambient_sim)
plotResiduals(sodium_amb_fec_no_atresia_lm)
testDispersion(sodium_amb_fec_no_atresia_lm)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.90609, p-value = 0.952
## alternative hypothesis: two.sided
# Normailty test
shapiro.test(residuals(sodium_amb_fec_no_atresia_lm))
##
## Shapiro-Wilk normality test
##
## data: residuals(sodium_amb_fec_no_atresia_lm)
## W = 0.86592, p-value = 0.08956
# Transform data: square root
sodium_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
mutate(Na...mmol.L. = sqrt(Na...mmol.L.))
# Rerun lm model with square root transformed data
sodium_amb_fec_no_atresia_lm_sqrt <- lm(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = sodium_ambient_sqrt)
summary(sodium_amb_fec_no_atresia_lm_sqrt)
##
## Call:
## lm(formula = Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = sodium_ambient_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.25662 -0.03796 0.02115 0.09537 0.12926
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.324e+01 8.923e-02 148.42
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 8.630e-05 5.083e-04 0.17
## Pr(>|t|)
## (Intercept) 6.31e-12 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.871
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1393 on 6 degrees of freedom
## Multiple R-squared: 0.004782, Adjusted R-squared: -0.1611
## F-statistic: 0.02883 on 1 and 6 DF, p-value: 0.8708
sodium_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = sodium_amb_fec_no_atresia_lm_sqrt)
plot(sodium_amb_fec_sqrt_sim)
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.75. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
plotResiduals(sodium_amb_fec_no_atresia_lm_sqrt)
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.75. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
testDispersion(sodium_amb_fec_no_atresia_lm_sqrt)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
sodium_amb_fec_no_atresia_lm_res <- residuals(sodium_amb_fec_no_atresia_lm_sqrt)
# Normality test
shapiro.test(sodium_amb_fec_no_atresia_lm_res)
##
## Shapiro-Wilk normality test
##
## data: sodium_amb_fec_no_atresia_lm_res
## W = 0.89593, p-value = 0.2654
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
#bartlett.test(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = sodium_ambient_sqrt)
# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = sodium_ambient_sqrt)
Sodium Ambient No Atresia Results
DHARMA Warning: Data dispersion error
Parametric assumptions: - Normality (Shapiro-Wilk test): Not Met - Variance (Bartlett’s test): NA
Nonparametric variance test (Levene’s test): NA
Residuals Plot (if significant)
# Na+
# Ambient samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
sodium_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Na...mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Sodium", x = "Weight Adjusted Fecudnity", y = "Sodium (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(sodium_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave(filename = "data-images/sodium_ambient_lm_scatterplot.pdf", plot = sodium_ambient_lm_scatterplot, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave(filename = "data-images/sodium_ambient_lm_scatterplot.png", plot = sodium_ambient_lm_scatterplot, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
Chloride Summary Stats
# Cl-
# Summary stats
# Ambient data
Ambient_Samples %>%
group_by(Ambient.Or.OAH, Consensus_Brood_Condition) %>%
summarize(count = n(),
median = round(median(Cl...mmol.L.), 3),
mean = round(mean(Cl...mmol.L.), 3),
sd = round(sd(Cl...mmol.L.), 3),
cv = round(sd(Cl...mmol.L.)/mean(Cl...mmol.L.), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 4 × 7
## Ambient.Or.OAH Consensus_Brood_Condition count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Ambient Excellent 1 158 158 NA NA
## 2 Ambient Normal 4 160 159. 3.86 0.024
## 3 Ambient Low 3 156 156. 0.577 0.004
## 4 Ambient Atresia 13 158 158. 4.61 0.029
Ambient_Samples %>%
group_by(Ambient.Or.OAH, Fecundity_Class) %>%
summarize(count = n(),
median = round(median(Cl...mmol.L.), 3),
mean = round(mean(Cl...mmol.L.), 3),
sd = round(sd(Cl...mmol.L.), 3),
cv = round(sd(Cl...mmol.L.)/mean(Cl...mmol.L.), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 3 × 7
## Ambient.Or.OAH Fecundity_Class count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Ambient High (>50,000) 4 158. 158 2.94 0.019
## 2 Ambient Low (~1,000s) 4 156 158. 3.70 0.023
## 3 Ambient Atresia 13 158 158. 4.61 0.029
Amb_Fec_No_Atresia_Samples %>%
group_by(Ambient.Or.OAH, Fecundity_Class) %>%
summarize(count = n(),
median = round(median(Cl...mmol.L.), 3),
mean = round(mean(Cl...mmol.L.), 3),
sd = round(sd(Cl...mmol.L.), 3),
cv = round(sd(Cl...mmol.L.)/mean(Cl...mmol.L.), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 2 × 7
## Ambient.Or.OAH Fecundity_Class count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Ambient High (>50,000) 4 158. 158 2.94 0.019
## 2 Ambient Low (~1,000s) 4 156 158. 3.70 0.023
Chloride Boxplot
# Cl- boxplot
# Ambient Samples
chloride_ambient_boxplot <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Consensus_Brood_Condition, y = Cl...mmol.L.)) +
geom_point(aes(x = Consensus_Brood_Condition, y = Cl...mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Chloride", x = "Parturition Success", y = "Chloride (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(chloride_ambient_boxplot)
ggsave(filename = "data-images/chloride_ambient_boxplot.pdf", plot = chloride_ambient_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/chloride_ambient_boxplot.png", plot = chloride_ambient_boxplot, device = "png")
## Saving 7 x 5 in image
chloride_ambient_fecundity_boxplot <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Cl...mmol.L.)) +
geom_point(aes(x = Fecundity_Class, y = Cl...mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Chloride", x = "Fecundity Class", y = "Chloride (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(chloride_ambient_fecundity_boxplot)
ggsave(filename = "data-images/chloride_ambient_fecundity_boxplot.pdf", plot = chloride_ambient_fecundity_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/chloride_ambient_fecundity_boxplot.png", plot = chloride_ambient_fecundity_boxplot, device = "png")
## Saving 7 x 5 in image
chloride_amb_fec_no_atresia_boxplot <- ggplot(data = Amb_Fec_No_Atresia_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Cl...mmol.L.)) +
geom_point(aes(x = Fecundity_Class, y = Cl...mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Chloride", x = "Fecundity Class", y = "Chloride (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(chloride_amb_fec_no_atresia_boxplot)
ggsave(filename = "data-images/chloride_amb_fec_no_atresia_boxplot.pdf", plot = chloride_amb_fec_no_atresia_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/chloride_amb_fec_no_atresia_boxplot.png", plot = chloride_amb_fec_no_atresia_boxplot, device = "png")
## Saving 7 x 5 in image
chloride_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
ggplot() +
geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Na...mmol.L.), colour = "black") +
labs(title = "Chloride", x = "Weight Adjusted Fecundity", y = "Chloride (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 20),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(chloride_amb_fec_no_atresia_scatterplot)
ggsave(filename = "data-images/chloride_amb_fec_no_atresia_scatterplot.pdf", plot = chloride_amb_fec_no_atresia_scatterplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/chloride_amb_fec_no_atresia_scatterplot.pdf", plot = chloride_amb_fec_no_atresia_scatterplot, device = "png")
## Saving 7 x 5 in image
Differences: Between Brood Conditon
# Na+
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
sodium_aov_ambient_brood <- aov(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
summary(sodium_aov_ambient_brood)
## Df Sum Sq Mean Sq F value Pr(>F)
## Consensus_Brood_Condition 3 176.8 58.95 0.622 0.611
## Residuals 17 1611.7 94.81
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(sodium_aov_ambient_brood)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## $Consensus_Brood_Condition
## diff lwr upr p adj
## Normal-Excellent 0.250000 -30.694633 31.19463 0.9999955
## Low-Excellent -3.666667 -35.626146 28.29281 0.9875929
## Atresia-Excellent 4.230769 -24.491760 32.95330 0.9745138
## Low-Normal -3.916667 -25.055875 17.22254 0.9514394
## Atresia-Normal 3.980769 -11.844573 19.80611 0.8897586
## Atresia-Low 7.897436 -9.830494 25.62537 0.5954090
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$Na...mmol.L.)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$Na...mmol.L.
## W = 0.52676, p-value = 3.625e-07
# QQplot:
sodium_ambient_res_qqplot <- ggqqplot(residuals(aov(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))) +
labs(title = "Sodium Interactive Residual QQplot",
subtitle = "residuals(aov(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))",
x = "Sodium Theoretical Quantiles (Predicted values)",
y = "Sodium Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(sodium_ambient_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
# bartlett.test(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
# Nonparametric variance test (Levene's test):
leveneTest(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 0.194 0.899
## 17
# Nonparametric stat test: Dunns
DunnTest(Na...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## Normal-Excellent 3.0000000 1.0000
## Low-Excellent -5.5000000 1.0000
## Atresia-Excellent 3.5769231 1.0000
## Low-Normal -8.5000000 0.3481
## Atresia-Normal 0.5769231 1.0000
## Atresia-Low 9.0769231 0.1252
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sodium Ambient Results
Parametric assumptions: - Normality (Shapiro-Wilk test): Not met - Variance (Bartletts test): NA
Nonparametric variance (Levene’s test): Met
Nonparametric Post Hoc: Dunn’s Test
Atresia-Low: mean.rank.diff = 9.0769231, pval = 0.1252
Differences: Between Fecundity Class
# Na+
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
sodium_aov_ambient_fec <- aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
summary(sodium_aov_ambient_fec)
## Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class 2 156.8 78.38 0.865 0.438
## Residuals 18 1631.8 90.66
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(sodium_aov_ambient_fec)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## $Fecundity_Class
## diff lwr upr p adj
## Low (~1,000s)-High (>50,000) -2.000000 -19.18271 15.18271 0.9526465
## Atresia-High (>50,000) 4.480769 -9.41330 18.37484 0.6939179
## Atresia-Low (~1,000s) 6.480769 -7.41330 20.37484 0.4737484
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$Na...mmol.L.)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$Na...mmol.L.
## W = 0.52676, p-value = 3.625e-07
# QQplot:
sodium_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Sodium Interactive Residual QQplot",
subtitle = "residuals(aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))",
x = "Sodium Theoretical Quantiles (Predicted values)",
y = "Sodium Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(sodium_ambient_fec_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
#bartlett.test(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
# Nonparametric variance test (Levene's test):
leveneTest(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.2098 0.8127
## 18
# Nonparametric stat test: Dunns
DunnettTest(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Dunnett's test for comparing several treatments with a control :
## 95% family-wise confidence level
##
## $`High (>50,000)`
## diff lwr.ci upr.ci pval
## Low (~1,000s)-High (>50,000) -2.000000 -18.008265 14.00826 0.9321
## Atresia-High (>50,000) 4.480769 -8.463634 17.42517 0.6108
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sodium Ambient Results
Parametric assumptions: - Normality (Shapiro-Wilk test): Not met - Variance (Bartlett’s test): NA
Nonparametric variance test (Levene’s test): Met
Correlations: Parameter to Numerical Fecundity (without Atresia IDs)
# Na+
# Linear Model
sodium_amb_fec_no_atresia_lm <- lm(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples)
summary(sodium_amb_fec_no_atresia_lm)
##
## Call:
## lm(formula = Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = Fecundity_No_Atresia_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.1021 -0.3010 0.6449 1.9545 3.3374
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 175.27869 2.08086 84.234
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.00558 0.01091 0.511
## Pr(>|t|)
## (Intercept) 4.4e-13 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.623
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.293 on 8 degrees of freedom
## (18 observations deleted due to missingness)
## Multiple R-squared: 0.03166, Adjusted R-squared: -0.08939
## F-statistic: 0.2615 on 1 and 8 DF, p-value: 0.6229
sodium_ambient_sim <- simulateResiduals(fittedModel = sodium_amb_fec_no_atresia_lm)
plot(sodium_ambient_sim)
plotResiduals(sodium_amb_fec_no_atresia_lm)
testDispersion(sodium_amb_fec_no_atresia_lm)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.90609, p-value = 0.952
## alternative hypothesis: two.sided
# Transform data: square root
sodium_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
mutate(Na...mmol.L. = sqrt(Na...mmol.L.))
# Rerun lm model with square root transformed data
sodium_amb_fec_no_atresia_lm_sqrt <- lm(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = sodium_ambient_sqrt)
summary(sodium_amb_fec_no_atresia_lm_sqrt)
##
## Call:
## lm(formula = Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = sodium_ambient_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.25662 -0.03796 0.02115 0.09537 0.12926
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.324e+01 8.923e-02 148.42
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 8.630e-05 5.083e-04 0.17
## Pr(>|t|)
## (Intercept) 6.31e-12 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.871
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1393 on 6 degrees of freedom
## Multiple R-squared: 0.004782, Adjusted R-squared: -0.1611
## F-statistic: 0.02883 on 1 and 6 DF, p-value: 0.8708
sodium_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = sodium_amb_fec_no_atresia_lm_sqrt)
plot(sodium_amb_fec_sqrt_sim)
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.75. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
plotResiduals(sodium_amb_fec_no_atresia_lm_sqrt)
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.75. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
testDispersion(sodium_amb_fec_no_atresia_lm_sqrt)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
sodium_amb_fec_no_atresia_lm_res <- residuals(sodium_amb_fec_no_atresia_lm_sqrt)
# Normality test
shapiro.test(sodium_amb_fec_no_atresia_lm_res)
##
## Shapiro-Wilk normality test
##
## data: sodium_amb_fec_no_atresia_lm_res
## W = 0.89593, p-value = 0.2654
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
#bartlett.test(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = sodium_ambient_sqrt)
# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = sodium_ambient_sqrt)
Sodium Ambient No Atresia Results
DHARMA Warning: Data dispersion error
Parametric assumptions: - Normality (Shapiro-Wilk test): Met - Variance (Bartlett’s test): NA
Nonparametric variance test (Levene’s test): NA
Residuals Plot (if significant)
# Na+
# Ambient samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
sodium_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Na...mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Sodium", x = "Weight Adjusted Fecudnity", y = "Sodium (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(sodium_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave(filename = "data-images/sodium_ambient_lm_scatterplot.pdf", plot = sodium_ambient_lm_scatterplot, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave(filename = "data-images/sodium_ambient_lm_scatterplot.png", plot = sodium_ambient_lm_scatterplot, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
Differences: Between Brood Condition
# Cl-
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
chloride_aov_ambient_brood <- aov(Cl...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
summary(chloride_aov_ambient_brood)
## Df Sum Sq Mean Sq F value Pr(>F)
## Consensus_Brood_Condition 3 24.59 8.197 0.463 0.711
## Residuals 17 300.65 17.685
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(chloride_aov_ambient_brood)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Cl...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## $Consensus_Brood_Condition
## diff lwr upr p adj
## Normal-Excellent 1.2500000 -12.114992 14.614992 0.9931731
## Low-Excellent -2.3333333 -16.136638 11.469971 0.9623849
## Atresia-Excellent 0.4615385 -11.943726 12.866803 0.9995589
## Low-Normal -3.5833333 -12.713361 5.546694 0.6851430
## Atresia-Normal -0.7884615 -7.623430 6.046506 0.9873943
## Atresia-Low 2.7948718 -4.861824 10.451567 0.7302966
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$Cl...mmol.L.)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$Cl...mmol.L.
## W = 0.96502, p-value = 0.6224
shapiro.test(residuals(aov(Cl...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov(Cl...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))
## W = 0.9502, p-value = 0.3437
# QQplot:
chloride_ambient_res_qqplot <- ggqqplot(residuals(aov(Cl...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))) +
labs(title = "Chloride Residual QQplot",
subtitle = "residuals(aov(Cl...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))",
x = "Chloride Theoretical Quantiles (Predicted values)",
y = "Chloride Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(chloride_ambient_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
# bartlett.test(Cl...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
# Nonparametric variance test (Levene's test):
leveneTest(Cl...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 1.3517 0.291
## 17
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Cl...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Cl...mmol.L. by Consensus_Brood_Condition
## Kruskal-Wallis chi-squared = 2.3847, df = 3, p-value = 0.4965
# Nonparametric Post Hoc: Dunns
DunnTest(Cl...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## Normal-Excellent 2.1250000 1.0000
## Low-Excellent -4.8333333 1.0000
## Atresia-Excellent 0.4615385 1.0000
## Low-Normal -6.9583333 0.8405
## Atresia-Normal -1.6634615 1.0000
## Atresia-Low 5.2948718 0.9031
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Chloride Ambient Results
Parametric assumptions: - Normality (Shapiro-Wilk test): Met - Variance (Bartletts test): NA
Nonparametric variance (Levene’s test): Met
Nonparametric Post Hoc: Dunn’s Test
Differences: Between Fecundity Class
# Cl-
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
chloride_aov_ambient_fec <- aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
summary(chloride_aov_ambient_fec)
## Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class 2 3.0 1.504 0.084 0.92
## Residuals 18 322.2 17.902
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(chloride_aov_ambient_fec)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## $Fecundity_Class
## diff lwr upr p adj
## Low (~1,000s)-High (>50,000) -0.5000000 -8.135556 7.135556 0.9847331
## Atresia-High (>50,000) 0.4615385 -5.712630 6.635707 0.9801564
## Atresia-Low (~1,000s) 0.9615385 -5.212630 7.135707 0.9170009
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$Cl...mmol.L.)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$Cl...mmol.L.
## W = 0.96502, p-value = 0.6224
shapiro.test(residuals(aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.96555, p-value = 0.6339
# QQplot:
chloride_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Chloride Residual QQplot",
subtitle = "residuals(aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))",
x = "Chloride Theoretical Quantiles (Predicted values)",
y = "Chloride Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(chloride_ambient_fec_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Cl...mmol.L. by Fecundity_Class
## Bartlett's K-squared = 0.82247, df = 2, p-value = 0.6628
# Nonparametric variance test (Levene's test):
leveneTest(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.5546 0.5838
## 18
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Cl...mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 0.4005, df = 2, p-value = 0.8185
# Nonparametric Post Hoc: Dunns
DunnettTest(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Dunnett's test for comparing several treatments with a control :
## 95% family-wise confidence level
##
## $`High (>50,000)`
## diff lwr.ci upr.ci pval
## Low (~1,000s)-High (>50,000) -0.5000000 -7.613663 6.613663 0.9777
## Atresia-High (>50,000) 0.4615385 -5.290623 6.213700 0.9711
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Sodium Ambient Results
Parametric assumptions: - Normality (Shapiro-Wilk test): Met - Variance (Bartlett’s test): Met
Nonparametric variance test (Levene’s test): Met
Correlations: Parameter to Numerical Fecundity (without Atresia IDs)
# Cl-
# Linear Model
chloride_amb_fec_no_atresia_lm <- lm(Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples)
summary(chloride_amb_fec_no_atresia_lm)
##
## Call:
## lm(formula = Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = Fecundity_No_Atresia_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.1972 -2.1453 -0.3044 2.3202 4.8677
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.581e+02 2.042e+00 77.415
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 7.028e-04 1.071e-02 0.066
## Pr(>|t|)
## (Intercept) 8.64e-13 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.949
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.232 on 8 degrees of freedom
## (18 observations deleted due to missingness)
## Multiple R-squared: 0.0005381, Adjusted R-squared: -0.1244
## F-statistic: 0.004307 on 1 and 8 DF, p-value: 0.9493
chloride_ambient_sim <- simulateResiduals(fittedModel = chloride_amb_fec_no_atresia_lm)
plot(chloride_ambient_sim)
plotResiduals(chloride_amb_fec_no_atresia_lm)
testDispersion(chloride_amb_fec_no_atresia_lm)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.90609, p-value = 0.952
## alternative hypothesis: two.sided
# Shapiro Test
shapiro.test(residuals(chloride_amb_fec_no_atresia_lm))
##
## Shapiro-Wilk normality test
##
## data: residuals(chloride_amb_fec_no_atresia_lm)
## W = 0.95118, p-value = 0.6824
# Transform data: square root
chloride_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
mutate(Cl...mmol.L. = sqrt(Cl...mmol.L.))
# Rerun lm model with square root transformed data
chloride_amb_fec_no_atresia_lm_sqrt <- lm(Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = chloride_ambient_sqrt)
summary(chloride_amb_fec_no_atresia_lm_sqrt)
##
## Call:
## lm(formula = Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = chloride_ambient_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.14856 -0.08721 -0.02279 0.07258 0.20200
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.257e+01 8.518e-02 147.566
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight -7.493e-05 4.852e-04 -0.154
## Pr(>|t|)
## (Intercept) 6.53e-12 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.882
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.133 on 6 degrees of freedom
## Multiple R-squared: 0.003959, Adjusted R-squared: -0.162
## F-statistic: 0.02385 on 1 and 6 DF, p-value: 0.8823
chloride_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = chloride_amb_fec_no_atresia_lm_sqrt)
plot(chloride_amb_fec_sqrt_sim)
## qu = 0.75, log(sigma) = -2.939317 : outer Newton did not converge fully.
plotResiduals(chloride_amb_fec_no_atresia_lm_sqrt)
## qu = 0.75, log(sigma) = -2.939317 : outer Newton did not converge fully.
testDispersion(chloride_amb_fec_no_atresia_lm_sqrt)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
chloride_amb_fec_no_atresia_lm_res <- residuals(chloride_amb_fec_no_atresia_lm_sqrt)
# Normality test
shapiro.test(chloride_amb_fec_no_atresia_lm_res)
##
## Shapiro-Wilk normality test
##
## data: chloride_amb_fec_no_atresia_lm_res
## W = 0.94043, p-value = 0.6153
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
#bartlett.test(Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = chloride_ambient_sqrt)
# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = chloride_ambient_sqrt)
Chloride Ambient No Atresia Results
DHARMA Warning: Data dispersion error
Parametric assumptions: - Normality (Shapiro-Wilk test): Met - Variance (Bartlett’s test): NA
Nonparametric variance test (Levene’s test): NA
Residuals Plot (if significant)
# Cl-
# Ambient samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
chloride_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Cl...mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Chloride", x = "Weight Adjusted Fecudnity", y = "Chloride (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(chloride_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave(filename = "data-images/chloride_ambient_lm_scatterplot.pdf", plot = chloride_ambient_lm_scatterplot, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave(filename = "data-images/chloride_ambient_lm_scatterplot.png", plot = chloride_ambient_lm_scatterplot, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
Potassium Summary Stats
# K+
# Summary stats
# Ambient data
Ambient_Samples %>%
group_by(Ambient.Or.OAH, Consensus_Brood_Condition) %>%
summarize(count = n(),
median = round(median(K...mmol.L.), 3),
mean = round(mean(K...mmol.L.), 3),
sd = round(sd(K...mmol.L.), 3),
cv = round(sd(K...mmol.L.)/mean(K...mmol.L.), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 4 × 7
## Ambient.Or.OAH Consensus_Brood_Condition count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Ambient Excellent 1 2.7 2.7 NA NA
## 2 Ambient Normal 4 2.7 2.7 0.082 0.03
## 3 Ambient Low 3 2.7 2.73 0.252 0.092
## 4 Ambient Atresia 13 2.7 2.86 0.304 0.106
Ambient_Samples %>%
group_by(Ambient.Or.OAH, Fecundity_Class) %>%
summarize(count = n(),
median = round(median(K...mmol.L.), 3),
mean = round(mean(K...mmol.L.), 3),
sd = round(sd(K...mmol.L.), 3),
cv = round(sd(K...mmol.L.)/mean(K...mmol.L.), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 3 × 7
## Ambient.Or.OAH Fecundity_Class count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Ambient High (>50,000) 4 2.7 2.7 0.082 0.03
## 2 Ambient Low (~1,000s) 4 2.7 2.72 0.206 0.076
## 3 Ambient Atresia 13 2.7 2.86 0.304 0.106
Amb_Fec_No_Atresia_Samples %>%
group_by(Ambient.Or.OAH, Fecundity_Class) %>%
summarize(count = n(),
median = round(median(K...mmol.L.), 3),
mean = round(mean(K...mmol.L.), 3),
sd = round(sd(K...mmol.L.), 3),
cv = round(sd(K...mmol.L.)/mean(K...mmol.L.), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 2 × 7
## Ambient.Or.OAH Fecundity_Class count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Ambient High (>50,000) 4 2.7 2.7 0.082 0.03
## 2 Ambient Low (~1,000s) 4 2.7 2.72 0.206 0.076
Potassium Boxplot
# K+ boxplot
# Ambient Samples
potassium_ambient_boxplot <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Consensus_Brood_Condition, y = K...mmol.L.)) +
geom_point(aes(x = Consensus_Brood_Condition, y = K...mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Potassium", x = "Parturition Success", y = "Potassium (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(potassium_ambient_boxplot)
ggsave(filename = "data-images/potassium_ambient_boxplot.pdf", plot = potassium_ambient_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/potassium_ambient_boxplot.png", plot = potassium_ambient_boxplot, device = "png")
## Saving 7 x 5 in image
potassium_ambient_fecundity_boxplot <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = K...mmol.L.)) +
geom_point(aes(x = Fecundity_Class, y = K...mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Potassium", x = "Fecundity Class", y = "Potassium (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(potassium_ambient_fecundity_boxplot)
ggsave(filename = "data-images/potassium_ambient_fecundity_boxplot.pdf", plot = potassium_ambient_fecundity_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/potassium_ambient_fecundity_boxplot.png", plot = potassium_ambient_fecundity_boxplot, device = "png")
## Saving 7 x 5 in image
potassium_amb_fec_no_atresia_boxplot <- ggplot(data = Amb_Fec_No_Atresia_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = K...mmol.L.)) +
geom_point(aes(x = Fecundity_Class, y = K...mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Potassium", x = "Fecundity Class", y = "Potassium (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(potassium_amb_fec_no_atresia_boxplot)
ggsave(filename = "data-images/potassium_amb_fec_no_atresia_boxplot.pdf", plot = potassium_amb_fec_no_atresia_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/potassium_amb_fec_no_atresia_boxplot.png", plot = potassium_amb_fec_no_atresia_boxplot, device = "png")
## Saving 7 x 5 in image
potassium_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
ggplot() +
geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = K...mmol.L.), colour = "black") +
labs(title = "Potassium", x = "Weight Adjusted Fecundity", y = "Potassium (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 20),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(potassium_amb_fec_no_atresia_scatterplot)
ggsave(filename = "data-images/potassium_amb_fec_no_atresia_scatterplot.pdf", plot = potassium_amb_fec_no_atresia_scatterplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/potassium_amb_fec_no_atresia_scatterplot.pdf", plot = potassium_amb_fec_no_atresia_scatterplot, device = "png")
## Saving 7 x 5 in image
Differences: Between Brood Condition
# K+
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
potassium_aov_ambient_brood <- aov(K...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
summary(potassium_aov_ambient_brood)
## Df Sum Sq Mean Sq F value Pr(>F)
## Consensus_Brood_Condition 3 0.1121 0.03736 0.505 0.684
## Residuals 17 1.2574 0.07397
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(potassium_aov_ambient_brood)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = K...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## $Consensus_Brood_Condition
## diff lwr upr p adj
## Normal-Excellent 0.00000000 -0.8643366 0.8643366 1.0000000
## Low-Excellent 0.03333333 -0.8593496 0.9260163 0.9995541
## Atresia-Excellent 0.16153846 -0.6407309 0.9638078 0.9389808
## Low-Normal 0.03333333 -0.5571209 0.6237876 0.9984688
## Atresia-Normal 0.16153846 -0.2804904 0.6035674 0.7296107
## Atresia-Low 0.12820513 -0.3669663 0.6233765 0.8812696
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$K...mmol.L.)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$K...mmol.L.
## W = 0.81347, p-value = 0.001062
shapiro.test(residuals(aov(K...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov(K...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))
## W = 0.886, p-value = 0.0189
# QQplot:
potassium_ambient_res_qqplot <- ggqqplot(residuals(aov(K...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))) +
labs(title = "Potassium Residual QQplot",
subtitle = "residuals(aov(K...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))",
x = "Potassium Theoretical Quantiles (Predicted values)",
y = "Potassium Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(potassium_ambient_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
# bartlett.test(K...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
# Nonparametric variance test (Levene's test):
leveneTest(K...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 0.6208 0.6111
## 17
# Nonparametric stat test: Kruskall Wallis
kruskal.test(K...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: K...mmol.L. by Consensus_Brood_Condition
## Kruskal-Wallis chi-squared = 1.2134, df = 3, p-value = 0.7498
# Nonparametric Post Hoc: Dunn's test
DunnTest(K...mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## Normal-Excellent 0.000000 1.0000
## Low-Excellent 0.500000 1.0000
## Atresia-Excellent 3.115385 1.0000
## Low-Normal 0.500000 1.0000
## Atresia-Normal 3.115385 1.0000
## Atresia-Low 2.615385 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Potassium Ambient Results
Parametric assumptions: - Normality (Shapiro-Wilk test): Not Met - Variance (Bartletts test): NA
Nonparametric variance (Levene’s test): Met
Nonparametric Post Hoc: Dunn’s Test
Differences: Between Fecundity Class
# K+
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
potassium_aov_ambient_fec <- aov(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
summary(potassium_aov_ambient_fec)
## Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class 2 0.1113 0.05563 0.796 0.466
## Residuals 18 1.2583 0.06990
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(potassium_aov_ambient_fec)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## $Fecundity_Class
## diff lwr upr p adj
## Low (~1,000s)-High (>50,000) 0.0250000 -0.4521380 0.5021380 0.9901955
## Atresia-High (>50,000) 0.1615385 -0.2242789 0.5473558 0.5449631
## Atresia-Low (~1,000s) 0.1365385 -0.2492789 0.5223558 0.6452587
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$K...mmol.L.)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$K...mmol.L.
## W = 0.81347, p-value = 0.001062
shapiro.test(residuals(aov(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.88368, p-value = 0.0171
# QQplot:
potassium_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Potassium Residual QQplot",
subtitle = "residuals(aov(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))",
x = "Potassium Theoretical Quantiles (Predicted values)",
y = "Potassium Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(potassium_ambient_fec_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: K...mmol.L. by Fecundity_Class
## Bartlett's K-squared = 4.6368, df = 2, p-value = 0.09843
# Nonparametric variance test (Levene's test):
leveneTest(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.7708 0.4773
## 18
# Nonparametric stat test: Kruskall Wallis
kruskal.test(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: K...mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 1.2081, df = 2, p-value = 0.5466
# Nonparametric stat test: Kruskall Wallis
kruskal.test(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: K...mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 1.2081, df = 2, p-value = 0.5466
# Nonparametric Post Hoc: Dunns test
DunnettTest(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Dunnett's test for comparing several treatments with a control :
## 95% family-wise confidence level
##
## $`High (>50,000)`
## diff lwr.ci upr.ci pval
## Low (~1,000s)-High (>50,000) 0.0250000 -0.4195254 0.4695254 0.9857
## Atresia-High (>50,000) 0.1615385 -0.1979081 0.5209850 0.4550
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Potassium Ambient Results
Parametric assumptions: - Normality (Shapiro-Wilk test): Not Met - Variance (Bartlett’s test): Kinda Met
Nonparametric variance test (Levene’s test): Met
Correlations: Parameter to Numerical Fecundity (without Atresia IDs)
# K+
# Linear Model
potassium_amb_fec_no_atresia_lm <- lm(K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples)
summary(potassium_amb_fec_no_atresia_lm)
##
## Call:
## lm(formula = K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = Fecundity_No_Atresia_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.21166 -0.11568 -0.04065 0.00147 0.39801
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 2.6473078 0.1245207 21.260
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.0007432 0.0006530 1.138
## Pr(>|t|)
## (Intercept) 2.52e-08 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.288
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1971 on 8 degrees of freedom
## (18 observations deleted due to missingness)
## Multiple R-squared: 0.1393, Adjusted R-squared: 0.03177
## F-statistic: 1.295 on 1 and 8 DF, p-value: 0.288
potassium_ambient_sim <- simulateResiduals(fittedModel = potassium_amb_fec_no_atresia_lm)
plot(potassium_ambient_sim)
## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully
plotResiduals(potassium_amb_fec_no_atresia_lm)
## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully
testDispersion(potassium_amb_fec_no_atresia_lm)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.90609, p-value = 0.952
## alternative hypothesis: two.sided
# Shapiro Test
shapiro.test(residuals(potassium_amb_fec_no_atresia_lm))
##
## Shapiro-Wilk normality test
##
## data: residuals(potassium_amb_fec_no_atresia_lm)
## W = 0.86151, p-value = 0.0795
# Transform data: square root
potassium_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
mutate(K...mmol.L. = sqrt(K...mmol.L.))
# Rerun lm model with square root transformed data
potassium_amb_fec_no_atresia_lm_sqrt <- lm(K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = potassium_ambient_sqrt)
summary(potassium_amb_fec_no_atresia_lm_sqrt)
##
## Call:
## lm(formula = K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = potassium_ambient_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.043790 -0.030608 -0.001939 0.010984 0.085420
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.6239747 0.0283420 57.299
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.0001535 0.0001614 0.951
## Pr(>|t|)
## (Intercept) 1.9e-09 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.378
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04424 on 6 degrees of freedom
## Multiple R-squared: 0.131, Adjusted R-squared: -0.01382
## F-statistic: 0.9046 on 1 and 6 DF, p-value: 0.3783
potassium_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = potassium_amb_fec_no_atresia_lm_sqrt)
plot(potassium_amb_fec_sqrt_sim)
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
## qu = 0.75, log(sigma) = -3.903504 : outer Newton did not converge fully.
plotResiduals(potassium_amb_fec_no_atresia_lm_sqrt)
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
## qu = 0.75, log(sigma) = -3.903504 : outer Newton did not converge fully.
testDispersion(potassium_amb_fec_no_atresia_lm_sqrt)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
potassium_amb_fec_no_atresia_lm_res <- residuals(potassium_amb_fec_no_atresia_lm_sqrt)
# Normality test
shapiro.test(potassium_amb_fec_no_atresia_lm_res)
##
## Shapiro-Wilk normality test
##
## data: potassium_amb_fec_no_atresia_lm_res
## W = 0.88746, p-value = 0.2216
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
# bartlett.test(K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = potassium_ambient_sqrt)
# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = potassium_ambient_sqrt)
Potassium Ambient No Atresia Results
DHARMA Warning: Data dispersion error
Parametric assumptions: - Normality (Shapiro-Wilk test): Kinda Met - Variance (Bartlett’s test): NA
Nonparametric variance test (Levene’s test): NA
Residuals Plot (if significant)
# K+
# Ambient samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
potassium_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = K...mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Potassium", x = "Weight Adjusted Fecudnity", y = "Potassium (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(potassium_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave(filename = "data-images/potassium_ambient_lm_scatterplot.pdf", plot = potassium_ambient_lm_scatterplot, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave(filename = "data-images/potassium_ambient_lm_scatterplot.png", plot = potassium_ambient_lm_scatterplot, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
Calcium Summary Stats
# Ca+2
# Summary stats
# Ambient data
Ambient_Samples %>%
group_by(Ambient.Or.OAH, Consensus_Brood_Condition) %>%
summarize(count = n(),
median = round(median(Ca....mmol.L.), 3),
mean = round(mean(Ca....mmol.L.), 3),
sd = round(sd(Ca....mmol.L.), 3),
cv = round(sd(Ca....mmol.L.)/mean(Ca....mmol.L.), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 4 × 7
## Ambient.Or.OAH Consensus_Brood_Condition count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Ambient Excellent 1 1.47 1.47 NA NA
## 2 Ambient Normal 4 1.42 1.43 0.05 0.035
## 3 Ambient Low 3 1.43 1.42 0.036 0.025
## 4 Ambient Atresia 13 1.24 1.27 0.114 0.09
Ambient_Samples %>%
group_by(Ambient.Or.OAH, Fecundity_Class) %>%
summarize(count = n(),
median = round(median(Ca....mmol.L.), 3),
mean = round(mean(Ca....mmol.L.), 3),
sd = round(sd(Ca....mmol.L.), 3),
cv = round(sd(Ca....mmol.L.)/mean(Ca....mmol.L.), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 3 × 7
## Ambient.Or.OAH Fecundity_Class count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Ambient High (>50,000) 4 1.45 1.45 0.044 0.03
## 2 Ambient Low (~1,000s) 4 1.41 1.41 0.033 0.023
## 3 Ambient Atresia 13 1.24 1.27 0.114 0.09
Amb_Fec_No_Atresia_Samples %>%
group_by(Ambient.Or.OAH, Fecundity_Class) %>%
summarize(count = n(),
median = round(median(Ca....mmol.L.), 3),
mean = round(mean(Ca....mmol.L.), 3),
sd = round(sd(Ca....mmol.L.), 3),
cv = round(sd(Ca....mmol.L.)/mean(Ca....mmol.L.), 3)) %>%
ungroup()
## `summarise()` has grouped output by 'Ambient.Or.OAH'. You can override using
## the `.groups` argument.
## # A tibble: 2 × 7
## Ambient.Or.OAH Fecundity_Class count median mean sd cv
## <ord> <ord> <int> <dbl> <dbl> <dbl> <dbl>
## 1 Ambient High (>50,000) 4 1.45 1.45 0.044 0.03
## 2 Ambient Low (~1,000s) 4 1.41 1.41 0.033 0.023
Calcium Boxplot
# Ca+2 boxplot
# Ambient Samples
calcium_ambient_boxplot <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Consensus_Brood_Condition, y = Ca....mmol.L.)) +
geom_point(aes(x = Consensus_Brood_Condition, y = Ca....mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Calcium", x = "Parturition Success", y = "Calcium (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(calcium_ambient_boxplot)
ggsave(filename = "data-images/calcium_ambient_boxplot.pdf", plot = calcium_ambient_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/calcium_ambient_boxplot.png", plot = calcium_ambient_boxplot, device = "png")
## Saving 7 x 5 in image
calcium_ambient_fecundity_boxplot <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Ca....mmol.L.)) +
geom_point(aes(x = Fecundity_Class, y = Ca....mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Calcium", x = "Fecundity Class", y = "Calcium (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(calcium_ambient_fecundity_boxplot)
ggsave(filename = "data-images/calcium_ambient_fecundity_boxplot.pdf", plot = calcium_ambient_fecundity_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/calcium_ambient_fecundity_boxplot.png", plot = calcium_ambient_fecundity_boxplot, device = "png")
## Saving 7 x 5 in image
calcium_amb_fec_no_atresia_boxplot <- ggplot(data = Amb_Fec_No_Atresia_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Ca....mmol.L.)) +
geom_point(aes(x = Fecundity_Class, y = Ca....mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Calcium", x = "Fecundity Class", y = "Calcium (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 18),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(calcium_amb_fec_no_atresia_boxplot)
ggsave(filename = "data-images/calcium_amb_fec_no_atresia_boxplot.pdf", plot = calcium_amb_fec_no_atresia_boxplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/calcium_amb_fec_no_atresia_boxplot.png", plot = calcium_amb_fec_no_atresia_boxplot, device = "png")
## Saving 7 x 5 in image
calcium_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
ggplot() +
geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Ca....mmol.L.), colour = "black") +
labs(title = "Calcium", x = "Weight Adjusted Fecundity", y = "Calcium (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 20),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(calcium_amb_fec_no_atresia_scatterplot)
ggsave(filename = "data-images/calcium_amb_fec_no_atresia_scatterplot.pdf", plot = calcium_amb_fec_no_atresia_scatterplot, device = "pdf")
## Saving 7 x 5 in image
ggsave(filename = "data-images/calcium_amb_fec_no_atresia_scatterplot.pdf", plot = calcium_amb_fec_no_atresia_scatterplot, device = "png")
## Saving 7 x 5 in image
Differences: Between Brood Condition
# Ca+2
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
calcium_aov_ambient_brood <- aov(Ca....mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
summary(calcium_aov_ambient_brood)
## Df Sum Sq Mean Sq F value Pr(>F)
## Consensus_Brood_Condition 3 0.1319 0.04396 4.527 0.0165 *
## Residuals 17 0.1651 0.00971
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(calcium_aov_ambient_brood)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Ca....mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## $Consensus_Brood_Condition
## diff lwr upr p adj
## Normal-Excellent -0.0400000 -0.3531867 0.2731867446 0.9830709
## Low-Excellent -0.0500000 -0.3734579 0.2734578789 0.9707679
## Atresia-Excellent -0.2007692 -0.4914663 0.0899278749 0.2401278
## Low-Normal -0.0100000 -0.2239473 0.2039472768 0.9991279
## Atresia-Normal -0.1607692 -0.3209355 -0.0006029264 0.0489650
## Atresia-Low -0.1507692 -0.3301914 0.0286529182 0.1173091
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$Ca....mmol.L.)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$Ca....mmol.L.
## W = 0.9463, p-value = 0.2894
shapiro.test(residuals(aov(Ca....mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov(Ca....mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))
## W = 0.91568, p-value = 0.07117
# QQplot:
calcium_ambient_res_qqplot <- ggqqplot(residuals(aov(Ca....mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))) +
labs(title = "Calcium Residual QQplot",
subtitle = "residuals(aov(Ca....mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples))",
x = "Calcium Theoretical Quantiles (Predicted values)",
y = "Calcium Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(calcium_ambient_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
# bartlett.test(Ca....mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
# Nonparametric variance test (Levene's test):
leveneTest(Ca....mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 0.9908 0.4207
## 17
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Ca....mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Ca....mmol.L. by Consensus_Brood_Condition
## Kruskal-Wallis chi-squared = 9.0798, df = 3, p-value = 0.02825
# Nonparametric Post Hoc: Dunn's test
DunnTest(Ca....mmol.L. ~ Consensus_Brood_Condition, data = Ambient_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## Normal-Excellent -3.125000 1.0000
## Low-Excellent -3.500000 1.0000
## Atresia-Excellent -11.153846 0.3321
## Low-Normal -0.375000 1.0000
## Atresia-Normal -8.028846 0.1412
## Atresia-Low -7.653846 0.2698
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Potassium Ambient Results
Which brood condition significantly different (if any)? - Atresia-Normal : (Tukey HSD: p adj = 0.0489650) - Atresia-Low : (Tueky HSD: p adj = 0.1173091)
Parametric assumptions: - Normality (Shapiro-Wilk test): Maybe Met - Variance (Bartletts test): NA
Nonparametric variance (Levene’s test): Met
Nonparametric Post Hoc: Dunn’s Test
Differences: Between Fecundity Class
# Ca+2
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
calcium_aov_ambient_fec <- aov(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
summary(calcium_aov_ambient_fec)
## Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class 2 0.1328 0.06641 7.281 0.00482 **
## Residuals 18 0.1642 0.00912
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(potassium_aov_ambient_fec)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## $Fecundity_Class
## diff lwr upr p adj
## Low (~1,000s)-High (>50,000) 0.0250000 -0.4521380 0.5021380 0.9901955
## Atresia-High (>50,000) 0.1615385 -0.2242789 0.5473558 0.5449631
## Atresia-Low (~1,000s) 0.1365385 -0.2492789 0.5223558 0.6452587
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$Ca....mmol.L.)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$Ca....mmol.L.
## W = 0.9463, p-value = 0.2894
shapiro.test(residuals(aov(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.91863, p-value = 0.08147
# QQplot:
calcium_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Calcium Residual QQplot",
subtitle = "residuals(aov(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples))",
x = "Calcium Theoretical Quantiles (Predicted values)",
y = "Calcium Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(calcium_ambient_fec_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Ca....mmol.L. by Fecundity_Class
## Bartlett's K-squared = 6.1296, df = 2, p-value = 0.04666
# Nonparametric variance test (Levene's test):
leveneTest(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 1.2704 0.3047
## 18
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Ca....mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 9.153, df = 2, p-value = 0.01029
# Nonparametric Post Hoc: Dunns test
DunnettTest(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Dunnett's test for comparing several treatments with a control :
## 95% family-wise confidence level
##
## $`High (>50,000)`
## diff lwr.ci upr.ci pval
## Low (~1,000s)-High (>50,000) -0.0375000 -0.1980658 0.1230658 0.7893
## Atresia-High (>50,000) -0.1807692 -0.3106040 -0.0509345 0.0071 **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Potassium Ambient Results
Parametric assumptions: - Normality (Shapiro-Wilk test): Maybe Met - Variance (Bartlett’s test): Not Met
Nonparametric variance test (Levene’s test): Met
Nonparametric Post Hoc: Kruskall Wallis test & Dunns test - Atresia-High (>50,000) : (Dunn’s test: pval = 0.0071 **)
Correlations: Parameter to Numerical Fecundity (without Atresia IDs)
# Ca+2
# Linear Model
calcium_amb_fec_no_atresia_lm <- lm(Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples)
summary(calcium_amb_fec_no_atresia_lm)
##
## Call:
## lm(formula = Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = Fecundity_No_Atresia_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.19339 -0.01787 0.01808 0.04987 0.10485
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.401e+00 5.882e-02 23.823
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight -3.752e-05 3.084e-04 -0.122
## Pr(>|t|)
## (Intercept) 1.03e-08 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.906
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09309 on 8 degrees of freedom
## (18 observations deleted due to missingness)
## Multiple R-squared: 0.001846, Adjusted R-squared: -0.1229
## F-statistic: 0.01479 on 1 and 8 DF, p-value: 0.9062
calcium_ambient_sim <- simulateResiduals(fittedModel = calcium_amb_fec_no_atresia_lm)
plot(calcium_ambient_sim)
## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully
plotResiduals(calcium_amb_fec_no_atresia_lm)
## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully
testDispersion(calcium_amb_fec_no_atresia_lm)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.90609, p-value = 0.952
## alternative hypothesis: two.sided
# Shapiro Test
shapiro.test(residuals(calcium_amb_fec_no_atresia_lm))
##
## Shapiro-Wilk normality test
##
## data: residuals(calcium_amb_fec_no_atresia_lm)
## W = 0.9084, p-value = 0.2702
# Transform data: square root
calcium_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
mutate(Ca....mmol.L. = sqrt(Ca....mmol.L.))
# Rerun lm model with square root transformed data
calcium_amb_fec_no_atresia_lm_sqrt <- lm(Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = calcium_ambient_sqrt)
summary(calcium_amb_fec_no_atresia_lm_sqrt)
##
## Call:
## lm(formula = Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = calcium_ambient_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.009302 -0.007894 -0.006232 0.004396 0.026830
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.180e+00 8.621e-03 136.816
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 1.143e-04 4.911e-05 2.328
## Pr(>|t|)
## (Intercept) 1.03e-11 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.0588 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01346 on 6 degrees of freedom
## Multiple R-squared: 0.4745, Adjusted R-squared: 0.3869
## F-statistic: 5.418 on 1 and 6 DF, p-value: 0.05882
calcium_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = calcium_amb_fec_no_atresia_lm_sqrt)
plot(calcium_amb_fec_sqrt_sim)
## qu = 0.25, log(sigma) = -3.52446 : outer Newton did not converge fully.
plotResiduals(calcium_amb_fec_no_atresia_lm_sqrt)
## qu = 0.25, log(sigma) = -3.52446 : outer Newton did not converge fully.
testDispersion(calcium_amb_fec_no_atresia_lm_sqrt)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
calcium_amb_fec_no_atresia_lm_res <- residuals(calcium_amb_fec_no_atresia_lm_sqrt)
# Normality test
shapiro.test(calcium_amb_fec_no_atresia_lm_res)
##
## Shapiro-Wilk normality test
##
## data: calcium_amb_fec_no_atresia_lm_res
## W = 0.77335, p-value = 0.01475
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
# bartlett.test(Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = calcium_ambient_sqrt)
# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = calcium_ambient_sqrt)
Potassium Ambient No Atresia Results
Parametric assumptions: - Normality (Shapiro-Wilk test): Met - Variance (Bartlett’s test): NA
Nonparametric variance test (Levene’s test): NA
Residuals Plot (if significant)
# Ca+2
# Ambient samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
calcium_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Ca....mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Calcium", x = "Weight Adjusted Fecudnity", y = "Calcium (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(calcium_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave(filename = "data-images/calcium_ambient_lm_scatterplot.pdf", plot = calcium_ambient_lm_scatterplot, device = "pdf")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
ggsave(filename = "data-images/calcium_ambient_lm_scatterplot.png", plot = calcium_ambient_lm_scatterplot, device = "png")
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
# Brood Condition Boxplots
ph_amb_bc <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Consensus_Brood_Condition, y = pH)) +
geom_point(aes(x = Consensus_Brood_Condition, y = pH), alpha = 0.5, colour = "black") +
labs(title = "pH", x = "Parturition Success", y = "pH") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
hct_amb_bc <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Consensus_Brood_Condition, y = Hct....)) +
geom_point(aes(x = Consensus_Brood_Condition, y = Hct....), alpha = 0.5, colour = "black") +
labs(title = "Hematocrit", x = "Parturition Success", y = "Hematocrit (%)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
glu_amb_bc <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Consensus_Brood_Condition, y = Glu..mmol.L.)) +
geom_point(aes(x = Consensus_Brood_Condition, y = Glu..mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Glucose", x = "Parturition Success", y = "Glucose (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
na_amb_bc <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Na...mmol.L.)) +
geom_point(aes(x = Fecundity_Class, y = Na...mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Sodium", x = "Fecundity Class", y = "Sodium (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
cl_amb_bc <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Consensus_Brood_Condition, y = Cl...mmol.L.)) +
geom_point(aes(x = Consensus_Brood_Condition, y = Cl...mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Chloride", x = "Parturition Success", y = "Chloride (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
k_amb_bc <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Consensus_Brood_Condition, y = K...mmol.L.)) +
geom_point(aes(x = Consensus_Brood_Condition, y = K...mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Potassium", x = "Parturition Success", y = "Potassium (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
ca_amb_bc <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Consensus_Brood_Condition, y = Ca....mmol.L.)) +
geom_point(aes(x = Consensus_Brood_Condition, y = Ca....mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Calcium", x = "Parturition Success", y = "Calcium (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
brood_condition_boxplots <- ph_amb_bc + hct_amb_bc + glu_amb_bc + na_amb_bc + cl_amb_bc + k_amb_bc + ca_amb_bc + plot_layout(ncol = 3, guides = "collect")
brood_condition_boxplots
# Fecundity class boxplots
ph_fec_class <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = pH)) +
geom_point(aes(x = Fecundity_Class, y = pH), alpha = 0.5, colour = "black") +
labs(title = "pH", x = "Fecundity Class", y = "pH") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
hct_fec_class <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Hct....)) +
geom_point(aes(x = Fecundity_Class, y = Hct....), alpha = 0.5, colour = "black") +
labs(title = "Hematocrit", x = "Fecundity Class", y = "Hematocrit (%)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
glu_fec_class <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Glu..mmol.L.)) +
geom_point(aes(x = Fecundity_Class, y = Glu..mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Glucose", x = "Fecundity Class", y = "Glucose (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
na_fec_class <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Na...mmol.L.)) +
geom_point(aes(x = Fecundity_Class, y = Na...mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Sodium", x = "Fecundity Class", y = "Sodium (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
cl_fec_class <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Cl...mmol.L.)) +
geom_point(aes(x = Fecundity_Class, y = Cl...mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Chloride", x = "Fecundity Class", y = "Chloride (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
k_fec_class <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = K...mmol.L.)) +
geom_point(aes(x = Fecundity_Class, y = K...mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Potassium", x = "Fecundity Class", y = "Potassium (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
ca_fec_class <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Ca....mmol.L.)) +
geom_point(aes(x = Fecundity_Class, y = Ca....mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Calcium", x = "Fecundity Class", y = "Calcium (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
fec_class_boxplots <- ph_fec_class + hct_fec_class + glu_fec_class + na_fec_class + cl_fec_class + k_fec_class + ca_fec_class + plot_layout(guides = "collect")
fec_class_boxplots
# lm plots
ph_am_lm <- ggplot(data = pH_ambient_sqrt, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = pH)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "pH", x = "Weight Adjusted Fecudnity", y = "Sqrt pH", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
hct_am_lm <- ggplot(data = hct_ambient_sqrt, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Hct....)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Hematocrit", x = "Weight Adjusted Fecudnity", y = "Sqrt Hematocrit (%)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
glucose_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Glu..mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Glucose", x = "Weight Adjusted Fecudnity", y = "Glucose (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
sodium_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Na...mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Sodium", x = "Weight Adjusted Fecudnity", y = "Sodium (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
chloride_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Cl...mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Chloride", x = "Weight Adjusted Fecudnity", y = "Chloride (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
potassium_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = K...mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Potassium", x = "Weight Adjusted Fecudnity", y = "Potassium (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
calcium_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Ca....mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Calcium", x = "Weight Adjusted Fecudnity", y = "Calcium (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
ambient_lm_plots <- ph_am_lm + hct_am_lm + glucose_am_lm + sodium_am_lm + chloride_am_lm + potassium_am_lm + calcium_am_lm + plot_layout(guides = "collect")
ambient_lm_plots
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).